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Abstract. We present a detailed analysis of the solution of the focusing nonlinear Schrodinger equation 
with initial condition ip{x,0) = N sech{x) in the limit N — » oo. We begin by presenting new and more 
accurate numerical reconstructions of the Af-soliton by inverse scattering (numerical linear algebra) for 
>rr^ ' N = 5, 10, 20, and 40. We then recast the inverse-scattering problem as a Riemann-Hilbert problem and 

' provide a rigorous asymptotic analysis of this problem in the large- Af limit. For those {x,t) where results 

have been obtained by other authors, we improve the error estimates from 0(N-'^^^) to 0{N-^). We also 
{^JQi analyze the Fourier power spectrum in this regime and relate the results to the optical phenomenon of 

' supercontinuum generation. We then study the Af-soliton for values of (x, t) where analysis has not been 

carried out before, and wc uncover new phenomena. The main discovery of this paper is the mathematical 
. mechanism for a secondary caustic (phase transition), which turns out to differ from the mechanism that 

generates the primary caustic. The mechanism for the generation of the secondary caustic depends essentially 
on the discrete nature of the spectrum for the Af-soliton, and more significantly, cannot be recovered from 
an analysis of an ostensibly similar Riemann-Hilbert problem in the conditions of which a certain formal 
, continuum limit is taken on an ad hoc basis. 
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1. Introduction 

It is well known (see, e.g., that the solution of the initial- value problem for the focusing nonlinear 
^ , Schrodinger equation 

g : (1) #t + + iv^iv = 

o 

(2) iP{x,0) = Nsech{x) 

is a special solution called the iV-soliton. This solution is rapidly decreasing in \x\ and periodic in t with 
period independent of N. To further describe this solution, we recall that (Q) is exactly solvable via the 
^ inverse-scattering framework introduced by Zakharov and Shabat in |16| . There are three steps in the 

procedure: 

(i) forward scattering — the initial data generates the scattering data which consists of eigenvalues, 
proportionality constants, and a reflection coefficient; 
■ (ii) time evolution — the scattering data have a simple evolution in time; 

(iii) inverse scattering — the solution of the partial differential equation at later times is reconstructed 
from the time-evolved scattering data. 

For the special initial data (0), Satsuma and Yajima [O] have shown that the reflection coefficient is iden- 
tically zero and there are N purely imaginary eigenvalues. Such a reflectionless solution whose eigenvalues 
have a common real part is an A^-soliton. 

Here, using the the inverse-scattering framework for (E^)-®, we study the A^-soliton in the limit N ^ oo. 
For the first few positive integer values of N, it is possible to write down explicit formulas for these exact 
solutions of the nonlinear partial differential equation iQ. However, these formulae rapidly become unwieldy. 
Even for = 3, the formula is already quite complicated and hard to analyze. Amazingly, while these 
formulae grow increasingly complicated as N increases, it is also true that certain orderly features emerge 
in the limit N —^ oo. 
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As a first step, we rescale -0 and t in CQl-lO to make tlie initial amplitude independent of TV and the 
period proportional to N, and we arrive at 

(3) ih^t + Y'P^^ + IV'I V = 0, 

(4) ipix,0) = Asech(x), 

where h = Hn '■= A/N and A> Q. Studying the large- limit of the A'-soliton is thus equivalent to studying 
the semiclassical [h J, 0) limit of (Pll^lB- A feature that emerges in the limit is the sharpening boundaries in 
the (x, t)-plane that separate different behaviors of the solution. Two such "phase transitions" were noticed 
in the numerical experiments of Miller and Kamvissis |12) . and the first, a so-called primary caustic curve in 
the (cc, <)-plane, was rigorously explained by Kamvissis, McLaughlin, and Miller |10j . Here, our main result 
is a description of the mechanism for the observed second phase transition; it differs from the first one. 

We use two formulations of the inverse-scattering step to study the limiting behavior of ©-O- When, as 
is the case here, the reflection coefficient is absent, the algebraic-integral system that one expects to solve to 
reconstruct the solution of for generic initial data reduces to a linear algebraic system. We derive such a 
system in §|21below. We also show, in S 13.11 that the reconstruction step can be cast as a discrete Riemann- 
Hilbert problem for a meromorphic 2x2 matrix unknown. We use the linear algebraic system as a basis for 
numerically computing the 7V-soliton while the Riemann-Hilbert formulation provides the starting point for 
our asymptotic analysis. In cither case, the reconstruction of the solution begins with the eigenvalues 

(5) Aw,fe (^fc+i^ ftw, k = Q,...,N-l, 
and the proportionality constants 

(6) 7JV,fc := k = Q,...,N-l. 

In § |31 we describe how to explicitly modify the discrete Riemann-Hilbert problem we obtain in § 13.11 
to arrive at an equivalent problem that is conducive to rigorous asymptotic analysis in the limit N ^ oo 
(or equivalently S — > 0). In S 13.21 we convert the discrete Riemann-Hilbert problem into a conventional 
Riemann-Hilbert problem for an unknown with specified jump discontinuities across contours in the complex 
plane. Next, we show how to "precondition" the resulting Riemann-Hilbert problem for the limit N oo 
by introducing a scalar function g{X) that is used to capture the most violent asymptotic behavior so that 
the residual may be analyzed rigorously. In some cases, the function g{X) leads to an asymptotic analysis 
in the limit N oo based on the theta functions of hypcrelliptic Riemann surfaces of even genus G, and 
in such cases the function g{X) can be constructed explicitly, as we show in S 13.41 This construction forces 
certain dynamics in {x, t) on the moduli (branch points) of this Riemann surface, and in S l3.4.1l we show that 
the moduli satisfy a universal system of quasilinear partial differential equations in Riemann-invariant form, 
the Whitham equations. The success or failure of the function g{\) constructed in this way in capturing the 
essential dynamics as A^ — > oo hinges upon certain inequalities and topological conditions on level curves that 
are described in S 13.4.21 With a function g{X) in hand that satisfies all of these conditions, we may proceed 
with the asymptotic analysis, which is based on a steepest descent technique for matrix Riemann-Hilbert 
problems developed by Deift and Zhou [7]. We summarize these steps in ^ 13.51 13.61 and 13.71 

The theory outlined above provides, through a handful of technical modifications, a refinement of results 
previously obtained by Kamvissis, McLaughlin, and Miller |1(J) . The key techincal improvements are the 
avoidance, through a dual interpolant approach developed in |11| and of a local parametrix near the 
origin in the complex plane and the explicit and careful tracking of the errors in replacing a distribution of 
point masses representing condensing soliton eigenvalues by its weak continuum limit, encoded in certain 
functions S{X) and T{X). These extra steps allow us to deduce the same asymptotic formulae obtained in 
|1U| governing the semiclassical limit up to and just beyond the primary caustic curve in the {x, i)-plane, 
but with an improved error estimate that is 0{h) (which beats the 0{h^^^) estimate in |1(J|). 

It is fair to say that much of the work in analyzing the A-soliton in the limit of large A^ comes about from 
vanquishing the poles representing the soliton eigenvalues from the matrix-valued Riemann-Hilbert problem 
that represents the inverse-scattering step. In some sense, the poles disappear with the interpolation step 
that converts the discrete Riemann-Hilbert problem with poles but no jump discontinuities into one with 
jump discontinuities and no poles. On the other hand, this is only a partial solution, since the jump matrix 
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relating the boundary values taken along a curve of jump discontinuity extends from this curve as an analytic 
function with poles at the soliton eigenvalues. However, as long as the jump contour remains bounded away 
from these poles, it is a reasonable approximation that can be controlled rigorously to replace the jump 
matrix by another one in which the distribution of poles is "condensed" into a continuous distribution with 
an analytic density. Making this replacement on an ad hoc basis changes the problem. We refer to this 
changed inverse-scattering problem as the formal continuum-limit problem and we discuss it briefly in § 13.81 
The formal continuum-limit problem indeed gives the correct asymptotics of the iV-soliton and related initial 
data for the semiclassical focusing nonlinear Schrodinger equation Q subject to the above caveat, and it 
even forms the starting point for the analysis of |15| . 

Interestingly, the choice of the function g{X) turns out to have the following two properties: 

• It depends only on the analytic weak limit of the discrete distribution of soliton eigenvalues. 

• It determines the contours of the Riemann-Hilbert problem of inverse scattering, and in particular 
their relation to the locus of accumulation of the discrete soliton eigenvalues. 

So, from the first property, g{X) docs not know about the poles of the jump matrix, and from the second 
property it chooses the contours appropriate for the asymptotic analysis. It is therefore possible that the 
contours selected by choice of g{X) could cross the locus of accumulation of discrete eigenvalues. And once 
the contours are no longer bounded away from the poles, the formal continuum-limit problem is not the 
correct model for inverse scattering. 

It turns out that, in the region beyond the primary caustic curve, this actually happens. In other words, 
the dependence of the function g{X) on {x, t) predicts a contour that as t increases passes through the branch 
cut that is the continuum limit of the pole distribution. At this point, the rigorous analysis must break 
down. Significantly, however, this is not the mechanism for further phase transitions. It turns out that 
the difficulty is a technical one that can be removed with the further use of the multi-interpolation method 
developed in and [Q. The additional steps that are required to surmount this crisis, and to therefore 
prove that there is no phase transition when curves determined from the function g{X) cross the pole locus, 
are described in 13.91 

The analysis described in § 13.91 is new, and one of the new features is the appearance of a new inequality 
that must be satisfied by the function g{X) on certain contours. It turns out that the failure of this inequality 
is the mathematical mechanism for the next phase transition of the iV-soliton, a secondary caustic curve in 
the {x, i)-plane. 

Before studying the secondary caustic, we discuss briefly in § 0] the quiescent region of the {x, t)-plane 
before the primary caustic curve. We recall in 14.11 some information from ^U] about the relation between 
the dynamics of the iV-soliton in this region and a family of hyperelliptic Riemann surfaces of genus zero. 
We also discuss the mathematical mechanism behind the primary caustic phase transition, which turns out 
to correspond to an instantaneous jump from genus zero to genus two. Then, in ^ 14. 21 we analyze the Fourier 
power spectrum of the A'^-soliton in the region before the primary caustic, and show that its evolution is 
consistent with the supercontinuum generation phenomenon of current interest in optical science (sec below). 

In carrying out a study of the A^-soliton for large A^, three different computational techniques come to 
mind: 

(i) Direct numerical simulation of the initial-value problem for the focusing nonlinear Schrodinger equa- 
tion ^ in the semiclassical limit. By adapting discretization methods to this problem, it is possible 
to study the dynamics of the A^-soliton as well as much more general initial data, and also non- 
integrable variants of (PJ. However, these methods are severely constrained in the limit of interest 
due to numerical stiffness, and worse yet, modulational instability that is exponentially strong in A'^. 
See mEllSlllI. 

(ii) Numerical solution of the inverse-scattering problem for the focusing nonlinear Schrodinger equation. 
This method applies only to the integrable equation Q, and like the direct numerical simulation 
method is also ill-conditioned in the semiclassical limit. However, it affords an important advantage 
over direct numerical simulation, namely that the solution is calculated independently for each (x, t), 
and therefore roundoff errors do not accumulate. 

(iii) Numerical construction of the function g{X). This method also applies only to the integrable equation 
Ipjjl. and it is further constrained in that it is only meaningful in the semiclassical limit. However, 
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it specifically takes advantage of the mathematical structure of this limit, and consequently the 
numerical calculation is extremely well-conditioned. 



In §0 we take the third approach and use numerical methods to solve for the function g(X) in the region 
between the primary and secondary caustics for the A''-soliton. While a numerical calculation, this is clearly 
fundamentally different from either the numerical solution of linear algebra problems equivalent to inverse 
scattering for the A^-soliton as in § or direct numerical simulation of the initial-value problem for (O 
corresponding to the iV-soliton. It is with the help of these calculations that we observe the crossing of the 
pole locus (which does not correspond to the secondary caustic, as can be seen once the analysis is modified 
as in S I3.9|I and ultimately the violation of the new inequality introduced in S 13. 91 f which causes the secondary 
caustic phase transition). 

In §|Hlwe take up the question of exactly what happens to the A'^-soliton immediately beyond the secondary 
caustic curve in the (a;, i)-plane. We characterize the secondary caustic (failure of the new inequality from 
S I3.9|I as the occurance of a critical point on a level curve defining the boundary of the region in the complex 
plane where the relevant inequality holds. We then make the guess that the failure of g{X) corresponding to a 
Riemann surface of genus two in this new way is resolved by going over to a formula for g{X) corresponding to a 
Riemann surface of genus four, but with significantly different topological features than in the (hypothetical) 
case that genus two fails due to the violation of the same inequality that leads to the primary caustic. We 
construct <?(A) in this new situation, and obtain an implicit description of the moduli of the corresponding 
genus four Riemann surface. These are new formulae that have a different character than those obtained 
earlier. We then prove that these new formula provide another solution to the same type of Whitham 
equations that govern the moduli in the case of simpler contour topology. 

We conclude the paper with some hypotheses regarding further phase transitions (higher-order caustics) 
for the A^-soliton and related initial- value problems for and by stressing once again that the phenomenon 
of the secondary caustic (and possibly further phase transitions) is fundamentally linked to the discrete nature 
of the eigenvalue distribution, which forces modifications to the analysis as described in S 13.91 and introduces 
coincident inequalities to be satisfied by g{X). These inequalities are simply not part of the asymptotic 
theory of the formal continuum-limit problem discussed in § 13.81 and therefore the latter problem is unable 
to correctly predict the secondary caustic. This came as some surprise to us, as we and other authors had 
always assumed that either the formal continuum-limit problem governs the dynamics for all time or that 
the secondary caustic occurs upon crossing the locus of eigenvalues. Both of these are incorrect. 

The problem Q for small h is relevant as a model for the propagation of light pulses in optical fibers that 
have the property of weak (because h is small) anomalous dispersion. There has always been great interest 
in optical fibers with weak dispersion. Initially this was because when such fibers are operated in the linear 
regime (small amplitude) the propagation is immune to dispersive spreading of pulses that is considered to 
degrade a data stream. For long fiber links, however, the theory of linear propagation becomes inadequate 
due to the accumulation of weakly nonlinear effects. More recently these fibers have become important 
again because it became clear that advantage could be taken specifically of the nonlinearity. Indeed, when 
operated in the weakly nonlinear regime (moderate amplitude) nonlinear effects are much stronger compared 
to linear effects and therefore propagation in such fibers can drastically alter the power spectrum of a signal, 
possibly in a useful way. Indeed, one of the applications envisioned for weakly dispersive fibers is so-called 
supercontinuum generation (see, for example, [HIIHI) in which a nearly monochromatic laser source is coupled 
into the fiber and spectrally broadened under propagation so that the output is a coherent source of white 
light. The output can then be filtered to produce coherent light of virtually any frequency, which is desirable 
for wavelength division multiplexing telecommunication systems. The most promising current technology 
for creating optical fibers with very low dispersion, in both the normal and anomalous regimes, is based 
on photonic crystal fibers, which are made from fiber preforms with extremely complicated cross sections. 
The complexity of the cross section is preserved to the microscopic level upon drawing the fiber (carefully 
stretching the cylindrical fiber preform along its axis until an optical fiber results), and as the possibilities for 
cross sections go far beyond the traditional corc/cladding/jacket step- index model it is possible to engineer 
fibers with properties thought impossible until recently [HI. 
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In this paper, we make frequent use of the three Pauh matrices: 



(7) 





■ 


1 ■ 




' -i' 




■ 1 " 


(Ti := 


1 





0-2 := 


i 


^3 := 


-1 



We use boldface notation {e.g. m or P) throughout for square matrices, with the exception of the identity 
matrix I, and arrow notation for column vectors {e.g. v, with row- vector transpose t/^). Complex conjugation 
is indicated with an asterisk, and P* means the matrix whose elements are the complex conjugates of the 
corresponding elements of P (no transpose). 



2. Linear Systems Describing the A^-Soliton 



2.1. Theory. The spectral data (0)-® can be used to set up a system of hnear equations whose solution 
yields the A^-soliton for arbitrary {x,t). (See, for example, the discussion in Chapter 2 of ^H].) To arrive at 
such a square inhomogeneous system, we begin by defining expressions 



(8) 



(9) 



and 



(10) 



N-l 

E 

p=0 



A{x,t,X) := V Ap{x,t)XP, 



C{x, t, A) + ^ Cp{x, t)XP, 

p=0 



Fjy^k{x, t) := exp ( - 2i(Ajv,fcX + X% ,^t)hN). 



Here, the unknown coefficients Aq{x, t), Ai{x,t), . . . , An^i{x, t) and Co{x,t),Ci{x,t), 
termined by the spectral data according to the relations 



, Cn-i{x, t) are de- 



(11) 
(12) 



A{x, t, XN,k)FN,k{x, t) = -fN,kC{x, t, XN,k), A: = . . . , TV - 1, 
C{x, t, X%^^)FN,k{x, ty = -j%kA{x, t, A^.,.), = 0, . . . , TV - 1. 



The A^-soliton solution of the focusing nonlinear Schrodinger equation is then given by 



(13) 



^N{x,t) ^ 2iAN-i{x,t). 



Clearly, the equations Hll|l and p2|l amount to a 2N x 2N system of linear equations for the coefficients 
Ap{x, t) and Cp{x, t), with coefficient matrices of block- Vandermonde type. One approach to computing the 
value of ipN{x,t) for fixed {x,t) is then to solve this linear algebra problem and then obtain ipiq{x^t) from 
(|13|l . This approach was used in U21 to find ipN{x,t) on a grid of values of {x,t) for = 5, = 10, and 

= 20. There is an advantage here over direct numerical simulation of the focusing nonlinear Schrodinger 
equation in that numerical errors do not propagate from one value of t to another. 

While the above approach is attractive, it can be made more so by reducing the problem explicitly to a 
linear algebra problem involving matrices of size only N x N. We use Lagrange interpolation and Hll|l to 
express C{x, t, A) in terms of the values A{x, t, XN,k), k = 0, . . . , N ~ 1. Thus, 



N-l 



(14) 



C{x,t,X) = X^ +J2 



Now we note that 
(15) 



Tl = 



N-l 



FN,n{x,t) 



A{x, t, XN,n) — X% „ 



n 



A- A 



N,j 



\N,n 



AjV, 



N-l 



A - A 



n=0 



J=0 



Xn.u ^ AtVj- 



JV-1 
J=0 



since the left-hand side is a monic polynomial of degree N which vanishes precisely at each the values 
Xn.o, ■ ■ ■ , Xn,n-i- Therefore, from (|14|l and H15(l we write 



N-l 



IN,', 



(16) C{x,t,X)^Yl 
We evaluate at and substitute into H12|l 



A{x,t,\N,n) 



N-1 

n 



A- A 



N-1 
J=0 



N-1 



(1^) E 



71=0 



lN,n 



-A{x,t,XN^n) 



N~l -V* A, 

j-j- ^ « i. ^ 



J=0 



N,k ^N.j 



N-1 
j=0 



-r^^,Aix,t,\*N,k). k = 0,...,N-l. 
Using the Lagrange interpolation formula once again, we write A in terms of its values at A^ o' ■ ■ ■ ' ^*n n-v 



N-l 



(18) 



A{x,t,X) = J2 



m=0 



N-1 



A{x,t,X% J J] 



A- A 



N.J 



\ * \ * 



Evaluating (|18|l at AAr,„ and substituting into H17() . we obtain 



N-l 



(19) 
where 



A{x, t, X^j.) + ^ WN,kmA{x, t, X^j.) = (j>N,k, k^Q,...,N ~l, 



ON,k 



PN,kj 



lN,k 

N-l 

-GN,k{x,t)* Yi i^*N,k - -^JVj); 
J=0 



N-l A, x* 



and 



N-l 



WN,km 



■= GN.k{x,t)* ^ GN,n{x,t)pN^ 



nmPN,kn- 



We denote by Wjv the N x N matrix with entries WN.km- The matrix Wat has a representation as 

Wat = B^Bjv, 

where Bj^jk ■= GN_j{x,t)p]\jjk- The linear system to solve is thus 

(20) {I + Wn)A = {1 + B%Bm)A^4, 

where A {A{x, t, X% q), . . . , A{x, t, A^ jy-i))^ ^-'^d ^ •= ('?^JV,o, ■ • ■ , 4'n,n-iV- Given a solution A of 
we recover the solution using via 

(21) ^N = 2iA-V, 
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where V := {Vna, ■ • ■ , Vn,n-i)'^ and 
(22) Vn. 

This follows from (fTH|) . since 

AT-l 

A(x,t,A) = ^ A(.T,t,A^,fc) 

N-l 

= ^ A(a:,t,A^^;,) 

2.2. Numerical implementation. In general, numerical integration of the initial-value problem 

for the semiclassically scaled focusing nonlinear Schrodinger equation is dificult; the problem is notoriously 
stiff due to multiple scales. The presence of small oscillations (microstructure) of wavelength and period 
of order h requires that one use a timestep proportional to h while resolving the larger scale structures 
(macrostructure) necessitates that the number of timesteps must be of order h^^. In addition, to accurately 
compute the spatial microstructure, the number of gridpoints/Fourier modes must also be proportional to 
h~^. Thus, simulating the semiclassical limit is computationally intensive, and the accumulation of roundoff 
errors is a serious issue. In spite of these difficulties, some careful numerical experiments have been carried 

out UliElE]. 

For the special case of initial data -0(2;, 0) = Asech(a;), the calculation described in the previous section 
provides a useful alternative approach. Rather than using some numerical integration scheme to directly 
compute an approximation to the solution of the partial differential equation, we may take advantage of 
the fact that the A^-soliton can be recovered from the solution of the linear algebraic system (|2UII . That 
is, for each fixed pair {x,t), we numerically solve 1)21) |l . An approach like this was first used by Miller and 
Kamvissis jT2] using instead the larger system of equations and (|12|) . 

In addition to being limited to particular initial data, the approach of solving (|20() suffers from the fact 
that the matrix I + W^r is ill-conditioned for large N and high-precision arithmetic is necessary for accurate 
computations. As the derivation leading to (|2()|l involves Lagrange interpolation on equi-spaced points, this 
is perhaps not surprising. On the other hand, an advantage of this approach is that the calculation for 
each (x, t) is independent of all the other calculations, so errors do not propagate in time and accumulate. 
Figures^^show density plots of \^n{x, <)p for TV = 5, 10, 20, and 40 computed by solving H20|l independently 
for a large number of (x, t) values using high-precision arithmetic. 

2.3. Phenomenology of the A^-soliton for large N. In each of Figures^^ the plotted solutions share 
some common features. After an initial period of smoothness, the solution changes over to a form with 
a central oscillatory region and quiescent tails. What is particularly striking in this sequence of figures is 
how the boundary between these two behaviors of the solution appears to become increasingly sharp as N 
increases. This initial boundary curve is called the primary caustic. The solution up to and just beyond the 
primary caustic has been carefully studied in |10j . See also the discussion in § 0] below. 

In Figure 121 and Figure 0] a second transition of solution behavior is clearly visible. It is this secondary 
caustic curve which is the main focus of this paper. As our subsequent analysis makes clear, this second phase 
transition is linked to the presence of the discrete soliton eigenvalues {AAr^fc} U {A^^.}, and the mechanism 
for this transition is different that of the transition across the primary caustic. 

We note that Tovbis, Venakides, and Zhou Jl] have studied the semiclassically scaled focusing nonlinear 
Schrodinger equation for a one-parameter family of special initial data. The forward-scattering procedure 
for this family 1jT4\ generates a reflection cocfficcnt and — depending on the value of the parameter — some 
discrete soliton eigenvalues. (By contrast, recall that our initial data is reflectionless.) For those values of the 
parameter for which there are no discrete eigenvalues, they prove that the solution undergoes only a single 
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N-l 

n 



A- A 



Na 



r. '^*N k~ ^N.;, 
n=0 ' 
n^k 



N-l 

n 



N-2\ 




0.0 0.5 1.0 1.5 2.0 

X 

Figure 1. The square modulus of ip^^x, t) plotted over a part of the positive quadrant of the 
{x, t) -plane. 

N=10 




X 



Figure 2. The square modulus of tpiQ^x^t) plotted over a part of the positive quadrant of 
the [x, t)-plane. 

phase transition and that there is no secondary caustic. When the value of the parameter is such that there 
are both soUtons and reflection, they leave the possibility of a second phase transition as an open question. 
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N=20 




X 




Figure 4. The square modulus of ipio^x^t) plotted over a part of the positive quadrant of 
the {x, t)-plane. 
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3. Method of Asymptotic Analysis 



3.1. Discrete Riemann-Hilbert problem. For each fixed positive integer N, and for fixed real values 
of X and t, consider solving the following problem: find a 2 x 2 matrix m(A; N, x, t) with entries that are 
rational functions of A such that 



• The poles are all simple and are confined to the points {^N.k}k=o ^^'^ {-^jv k}k=o^ such that 



(23) 



(24) 



Res m(A; A^, x, t) 



lim m(A; A^, x', t) 





CN,k{x,t) 



and 



Res m(A; Af, t) = lim m(A;A, a;,i) 



hold for fc = 0, 



-CAr.fc(x,t)* 





, A - 1, where CN,k{x,t) ■= c?^ ^,e2''3(^"-'='^'*)/'*«, and 



(25) 



^N,k 



N-1 

Ha 

n=0 



N,k 



lN,k 



N-1 

Ha 



Q{X;x,t) := Xx + X^t . 



N,k 



The matrix m(A; N, x, t) is normalized so that 



(26) 



lim m{X]N,x,t) = I. 



Then, the function defined in terms of m(A;A^, x,t) by the hmit 



(27) 



ipNix,t) := 2i lim Ami2(A; A, x, t) 

A — ^oo 



is the A-soliton solution of the semiclassically scaled focusing nonlinear Schrodinger equation. 

This Riemann-Hilbert problem essentially encodes the linear equations introduced in § |21 describing the 
A-soliton, in a way that is conducive to asymptotic analysis in the limit A — > cx). That the function 
tpN^x, t) solves the semiclassically scaled focusing nonlinear Schrodinger equation (O is easy to show. Indeed, 
let n(A;A, x,t) = m(A; A^, i)e~''5(-'';a:,t)o-3/?iN^ Then, it is easy to check that the residue conditions on 
m(A; N, x, t) translate into analogous conditions on n(A; A, .t, t) that are independent of x and t: 



(28) 
and 
(29) 



Res n(A; AT, x, t) = hm n(A; AT, x, t) 



Res n(A; AT, x, t) = hm n(A; A^, x, t) 







-Art 











It is easy to see that det(n(A; A^, x, i)) = 1, because clearly the determinant is a meromorphic function, 
with possible simple poles only at the points {XN,k} U {X*j^ j.}, that tends to one as A ^ oo. But from the 
residue conditions it is easy to verify that det(n(A; N, x, t)) is regular at the possible poles, and so is an entire 
function tending to one at infinity that by Liouville's Theorem must be constant. In particular, n(A; x, t) 
is nonsingular for all A, so it follows that dxri{X; N, x,t)n{X; N, x,t)~^ and dtn{X; N, x,t)n{X; N, x,t)~^ are 
both entire functions of A. Moreover, they are polynomials in A, as can be deduced from their growth at 
infinity. Indeed, m(A; A, x, t) necessarily has an expansion 



(30) 



m(A; A, X, i) = I + ^ A ^nip , nip — mp{N, x,t) , 
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which is uniformly convergent for |A| sufficiently large (larger than the modulus of any prescribed singularity 
of in{X] N,x,t) is enough), and differentiable term-by-term with respect to x and t. Therefore, 



d^n{X; N, x, t)n{X; N, x, 



dx 



-iQ{\:x,t)a3/hK 



(31) 



Similarly, 



(l + A ^mi + A + ■ ■ 
e^Qi^'^,t)'^^/'iN (J _ x-^mi + A-2(m2 - ma) + • • • ) 

[-ihJf^Xas - ifi^^miCTs H ) • (I - A^-^mi H ) 

-ihjj^ (Acts + [1111,(73]) . 



dtn{X;N,x,t)n{X;N,x,ty^ 



dt 



-iQ{\;x,t)a3/hN 



(32) 



(1 + A ^mi + X^'^m.2 + ■ ■ 
^^Q{\;x,t)a3/h« (I - X-^mi + X-'^{m\ - ma) + 

'-ih'^^X^cr's, - i^^^AmiCTs - i^i^^maCTa H ) 

(I-A-imi+A-2(m2-m2) + ---) 



= -ift^^ (AV3 + A[mi, (73] + [m2, a-i] - [mi, cr3]mi) . 

Consequently, n(A; iV, x, t) is a simultaneous fundamental solution matrix for general A of the linear differ- 
ential equations 

ihNdxn{X;N,x,t) = (Atrs + [mi, 0-3]) n(A; iV, x, t) 

(33) 

ihNdtn{X;N,x,t) = (A^cr3 -(- A[mi, 0-3] -I- [m2, 0-3] - [mi, CTsJmi) n(A; iV, x, t) . 
The coefficient matrices therefore satisfy the zero-curvature compatibility condition 

ifiNdt {Xa3 + [mi, (T3]) - ifiNdx {X^cr-i + A[mi, (73] + [m2, (73] - [mi, cr3]mi) 

(34) 

+ [Act3 + [mi,(73], AV3 + A[mi,(73] + [m2,(73] - [mi,(73]mi] = 0. 
Separating out the coefficients of the powers of A wc obtain two nontrivial equations: 

(35) -ihNdx[nii,a3] - [[m2, (73], (73] + [[mi, (73] mi, (73] ^ 0, 
and 

(36) ihNdt[nii,a3] - ih^dx ([m2,cr3] - [mi,(73]mi) + [[mi, (73], [m2,(73]] - [[mi, (73], [mi,cr3]mi] 
Introducing the notation 

(37) A = + A^^ , 

for separating a 2 x 2 matrix into its diagonal and off-diagonal parts, wc have for any A: 

(38) [A,a3] =2A0^(73, 
The first equation becomes 

(39) -iftAT^^mf ^(73 - 2m°^ + 2mf ^mf = , 
which is purely off-diagonal, while the second equation has both diagonal parts: 

(40) zft^9,(mf ^agmf ^) + 2[m?^a3, m0^a3] - 2[m?^a3, mf ^c73mf ] 
and off-diagonal parts: 

(41) ^hr,^,{rar^,) ~ '-^^dx{m<^^ a,) + 



0. 



a,(mf^mf (73) - 2[m^-a3,mr^3m^-] = 



on. 



0, 



Eliminating rrij using the first equation, the ofT-diagonal part of the second equation becomes: 



(42) 



4m' 



OD3 



In other words, if we introduce notation for the off-diagonal elements of mi as follows, 



(43) 



m 



OD 



q 
r 



then q and r satisfy a coupled system of partial differential equations: 

(44) ihNdtq+ ^dlq-Arq^ = -ifiNdtr + ^d^r ~ Aqr"^ = Q . 



If for all X and t we have 



(45) 



q*, then these become the focusing nonlinear Schrodinger equation 

2ir* . 



f{2 



Therefore, to complete the proof, it remains only to show that indeed r = — g*. To do this we consider 
along with the solution m(A; N, x, t) of the discrete Riemann-Hilbert problem the corresponding matrix 
m.{\; N ,x,t) := (T2m(A*; iV, x, t)*(T2, where the star denotes componentwise complex conjugation. Clearly 
m(A; x, t) is analytic in A with simple poles at {Ajv^fcjU {AJ^ j.} (because the pole set is complex-conjugate 
invariant) and tends to the identity as A ^ oo (because ct| = I)- Furthermore, 



Res m(A; N, x, t) 



a2 Res m(A; A^, X, t) a2 

\ i, / 



(46) 



(72 liin m{X;N,x,t) 



lim m(A; N, x, t)(j2 





-CN.k{x, t) 





0-2 



0-2 



= lim m(A; x, t) 





CN,k{x,t) 



By a similar calculation, 

(47) Res m(A; iV, x, t) = lim m(A; N, x, t) 



~CN,k{x,t)* 





It follows that m(A; A^, X, t)m(A; iV, X, t)~^ is an entire function of A that tends to the identity matrix as 
A oo. By Liouville's Theorem, we thus have m(A; A^, x,t) = m(A;iV, x,t). Expanding both sides of this 
identity near A = oo, we get 



(48) 

so in particular, 
(49) 



I + A-^mi = I + A" V2mtcr2 + 0(A" 



nil = cr2nii (72 ■ 



which gives r = —q*. 

That the function ipN{x,t) satisfies iPn{x,0) = Asech(x) for all N is more difficult to show by studying 
properties of the matrix m(A; N, x, t). In general, this follows from solving the corresponding direct-scattering 
problem which was done with the help of hypergeometric functions by Satsuma and Yajima [1 3 j . Here we 
illustrate the corresponding inverse-scattering calculation in the most tractable case of = 1. When A^ = 1, 
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m(A; 1, X, 0) may be sought in the form 
(50) m(A;l,a;,0) 



A + au{x) ai2{x) 
X-iA/2 

X-iA/2 



The residue relations then say that 

iA/2 + aii{x) = -iAe^"" 



(51) 



a2i{x) 
ai2ix) 



-iAe-"" 



-iAe 



X + iA/2 

A + a22{x) 
X + iA/2 

lA ' 

iA/2 + a22{x) 
iA 

iA/2 - aii(x) 
iA 

-a2i{x) 



iA 



~iA/2 + a22{x) = -iAe^"" 
It foUows that ai2{x) ~ — sech(.T)/2, from which we indeed find that ■!/'i(2;, 0) = ^sech(a;). 
3.2. First modification: removal of poles. 

3.2.1. Interpolation of residues. In the foUowing, to keep the notation as simple as possible, we suppress 
the parametric dependence on N, x, and t. We first modify the matrix unknown m(A) := m(A;iV, x,t) by 
multiplication on the right by an explicit matrix factor which differs from the identity matrix in the regions 
Di, -D-i, and their complex conjugates, as shown in Figure [S] Let ii'o denote the common boundary arc 




Figure 5. The regions Di and -D-i in the upper half-plane, and the oriented boundary arcs 
K^i, Kq, and Ki. The N soliton eigenvalues {XN,k}k=o '^^'^ imaginary interval [0, i^] 
in which they accumulate as N oo are shown in red. The number e > will be taken later 
to be sufficiently small, but independent of N. 



of Di and D-i oriented in the direction away from the origin. Let Ki denote the remaining part of the 
boundary of Di lying in the open upper half-plane, oriented in the direction toward A = e. Let denote 
the remaining part of the boundary of -D_i lying in the open upper half-plane, oriented in the direction 
toward A ~ — e. We will call the point where the three contour arcs Kj meet Aq. 
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Note that the region D^i contains the sohton eigenvalues {XN,k}k=o "^^^ region Di is needed 

for a technical reason; with its help we will be able to ultimately remove some jump discontinuities from the 
neighborhood of A = 0. We set 



(52) 



M(A) := m(A) 



1 

-iP(A)e[2'Q(^)+»«°(^)l/'i 1 



for X e Di, 



(53) 



M(A) := m(A) 



1 



for A e D^i, 



For all A in the upper half-plane outside the closure of Di U ZJ-i, we set ]V[(A) = m(A). For A in the lower 
half-plane, we define M(A) := (72M(A*)*(T2, where the star denotes componentwise complex-conjugation. 
Here, we are using the notation 

(54) P(A) := n 



fe=0 



A- A 



N,k 



and 
(55) 

and, h ~ h^. It is easy to see that the matrix M(A) is holomorphic at the soliton eigenvalues A — XN,k 
where m(A) has its poles in the upper half-plane. 



^(A) := inX + ttA, 



3.2.2. Aside: even symmetry of ipN{x,t) in x and the formal continuum limit. Since M(A) = m(A) in a 
neighborhood of A = oo, the iV-soliton is equivalently defined by the formula 



(56) 



2i lim AA/i2(A) 

A — ^oo 



From the conditions determining M(A) is is easy to show that V'7v(a;,t) is, for each N and each t, an even 
function of x. For this purpose, we may suppose without any modification of ipi\jix,t) that e = 0, that 
Di = 0, and that D^i is symmetric about the imaginary axis. We also temporarily re-introduce the explicit 
parametric dependence on x, and suppose that D^i is fixed as a; S M varies. Then the matrix M(— A; x) has 
the same domain of analyticity as M(A; x), and therefore we may compare M(A; x) with the matrix Ml* (A; x) 
defined by setting 

^^-[2iiXx-\^t)+ie°ix)]/h P{Xy^ 



(57) M«(A; x) cr3M(-A; x)a3 



X e 



while elsewhere in the upper half-plane we set 

(58) M«(A; x) := a^Mi-X; x)a3P{XY' , 

and then we define Ml'(A; x) for A in the lower half-plane by setting M'*(A; x) := f72M*(A*; x)*f72 for 3(A) < 0. 
Note that because the zeros of P{X) are confined to D'^i while the poles are confined to D-i, this defines 
M'(A; x) as a sectionally holomorphic function of A with the same (cc-independent) domain of analyticity as 
M(A;x). 

Now, if we use a subscript (respectively ") to denote a boundary value taken on the boundary of 
_D_i from inside (respectively outside), then it is an easy exercise to check that 

(59) Mi{X;xy^Ml{X;x) = M_(A; -a-)-^M+(A; -x) , for x e R and A £ dD^iUdDl-^ . 

Also, since P{X) ^ 1 as A ^ oo, it follows that M'*(A;x) ^ I as A ^ oo for each a; e R. By Liouville's 
Theorem, which assures the uniqueness of M(A;x) given the jump condition across dD^i U dD*_T^ and the 
normalization condition at A = cxj, it follows that 



(60) 

Using this relation, we clearly have that 



M''(A;a;) =M(A;-a;). 



(61) 



2i lim XMl^{X]x) 

X — 'oc 



2i lim AA/i2(A; — x) = 7/'Ar(— X, i) , 

A— i-DO 



where the last equality follows from H56|l . On the other hand, directly from the definition H58|) valid for 
sufficiently large |A|, we may use the fact that P{X) ^ 1 as A — > oo to conclude that 



(62) 



2i lim AAf»2(A;x-) 



-2i lim AMi2(-A;a;) 

A — yoo 



2i lim fiMi2{p-; x) — iPn{x, t) . 



Therefore we learn that iPn{—x, t) = iPn{x, t) holds for all real x, so the A^-soliton is an even function of x for 
each t. Therefore, in all remaining calculations in this paper we will suppose without loss of generality that 
a; > 0. (Since from we see that time reversal is equivalent to complex conjugation of '0jv(a^,f) because 
iPn{x,0) is real, we will also assume that t > 0.) These choices lead to certain asymmetries in the complex 
plane as already apparent in Figure |S1 for a discussion of how choice of signs of x and t relates to the parity 
of certain structures in the complex spectral plane that are introduced to aid in asymptotic analysis, see 

nni- 

The even symmetry of V'jv(a;,t) in x is not an earth-shattering result, of course, but what is interesting 
is that the same argument fails completely if a natural continuum limit related to the limit ^ oo is 
introduced on an ad hoc basis. Indeed, it is natural to consider replacing the product P{X) by a formal 
continuum limit by "condensing the poles" as follows: define 



(63) 



P(A) := cxp 



iA 



log(— i(A — 77)) rfr/ — i / log{—i{X — 1])) drj 



which may be viewed as coming from interpreting the sums in the exact formula 



(64) 



P(A) = cxp 



fe=0 



7V-1 

E 

k=Q 



log(A - XN,k)fi 



as Riemann sums and passing to the natural integral limit for A fixed. The function P{X) is analytic for 
A e C \ [—1^,1^], and where P(A) has accumulating poles and zeros, P(A) has a logarithmic branch cut. 
Then, while M(A; x) is defined as being a matrix with the symmetry M(A; x) = cr2M(A*; x)*cr2 that satisfies 
the normalization condition M(A; x) — > I as A — > 00 and is analytic except on dD^i U dD*_i along which it 
takes continuous boundary values related by 



(65) 



M+(A;a;) = M_(A;a;) 



1 

iP(A)e[2*(^^+^'*)-*''"(^)l/'^ 1 



A e dD^ 



we may define a matrix M(A) that satisfies exactly the same conditions as M(A) but with P(A) replaced by 
P(A) in the jump condition. It is a direct matter to show that like iPn{x, t), the function defined by 



(66) 



ipNiXjt) := 2i lim AA/i2(A;x) 



is also a solution of the semiclassically scaled focusing nonlinear Schrodinger equation (virtually the same 
arguments apply as was used to prove this about ipN{x^t)). However, whether il)iq(x,t) is an even function 
of x is in our opinion an open question. Indeed, if we try to mimic the above proof of evenness of iPn{x, t), 
we would be inclined to try to define a matrix M*'(A; x) starting from M(A; x) by formulae analogous to (|57() 
and (|58|l . but with P(A) replaced everywhere by P(A). Comparing M*(A; x) with M(A; —x) then becomes a 
problem, because while M(A; ~x) is analytic in D-i U -Dli, the matrix M.^[X;x) has a jump discontinuity 
across the segment [— iA,?^]. Thus, in stark contrast with H6()(l wc have 

(67) M''(A;a;) 7^M(A;-x), 

and we cannot conclude (by a completely analogous proof, anyway) any evenness of ip{x,t). 

The possibility that ip]\f{x, t) may not be an even function of x while ip]^{x, t) most certainly is even casts 
some doubt on the prospect that TpN(x,t) might be a good approximation to ipN(x,t). As it is TpNix,t) 
that is related to the solution of the "continuum- limit" Ricmann-Hilbcrt problem for a matrix P(A) to be 
introduced in § 13.81 we have some reason to suspect at this point that without careful accounting of the 
errors introduced by making the ad hoc substitution P(A) — > P(A), a study of the "condensed-pole" problem 
may not be relevant at all to the asymptotic analysis of the A'^-soliton '4>n{x, t), at least for certain x and t. 
We will give evidence in this paper that such suspicion is entirely justifiable. 
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3.3. Second modification: introduction of g-function. Suppose that g(A) is a function analytic for 
A e C \ (A'_i U A'o U Kli U A'o) tliat satisfies the symmetry condition 

(68) 5(A) + .g(A*)* = 0, 

and g{X) — > as A ^ oo. Note that in particular g{X) is analytic in the real intervals (— e, 0) and (0, e). We 
change variables again to a matrix function N(A) by setting 

(69) N(A) :=M(A)e-s(^)'"='/''. 

Letting the subscript "+" (respectively "— ") denote a boundary value taken on one of the contours Kj from 
the left (respectively right) according to its orientation, we deduce from our definitions and the continuity 
of m(A) across each of these contours that the following "jump relations" hold: 

iP(A)e[2*'3(A)-»e°(A)-g+(A)-s_(A)]/R g[s+(A)-g-(A)]/ft 



(70) 



N+(A) = N_(A) 



A G X_i , 



(71) 



N+(A) = N_(A) 



1 

iP(A)e[2'Q(^)+*«"(^)-25(A)]/'i 1 



A e i^i 



(72) N+(A)=N_(A) 



The jump relations holding on the conjugate contours in the lower half-plane follow from the symmetry 
N(A*) = cr2N(A)*(T2. Finally, there are also jump discontinuities across the real intervals (— e, 0) and (0,e), 
both of which we assign an orientation from left to right. Using the above symmetry relation along with the 
facts (holding for x and t real and h = A/N with N E "Z): 

(73) P{X*y = P(A)-i , 0(A*)* = Q(A) , 0"{X*)* = -0"{X) + 27ThN , 

we find that 
(74) 



N+(A) = N_(A) 



1 



-2ieo{\)/h 



iP{X)e 



[2iQ{\)-ie°(\)-2g(\)]/h 



'iP(A)~ie[~2iQ(A)-je«(A)+2s(A)]/ri 



(75) N+(A)=N_(A) 



1 + e2*''"(^)/'^ iP(A)- 
-iP(A)e[2«QW+»«"W-23(A)]/;i 



1 



-2iQ(\) + ie"{X)+2g{\)]/h 



Ae (-6,0) 



AG (0,e), 



(76) 



L+(A)-L_(A) = 



Let C denote a simple contour lying in D-i starting from the origin and terminating at A = lA. There is 
a unique function L{X) with the properties that (i) it is analytic for A e C\ (CUC*), (ii) it takes continuous 
boundary values on each side of the contour C U C* satisfying 

-2ie^{X) , for A e C, 

-2ie^{X*)* , for A e C\ 

where i+(A) (respectively L_(A)) refers to the boundary value taken on C U C* from the left (respectively 
right) as it is traversed from —iA to lA^ and (iii) it is normalized so that L{X) tends to zero as A ^ oo. 
Indeed, the function LiX) is easily seen to be unique from these conditions (by Liouvillc's theorem and 
continuity of the boundary values as is compatible with the conditions (|76|l ') and we will give an explicit 
construction later (see (|102|l - Hll)7|l '). The function L{X) enjoys the following symmetry property: 

(77) L(A) + P(A*)* =0. 

Related to L{X) is another function L{X) defined as follows. Let Coo denote an infinite simple contour in 
the upper half-plane emanating from X ~ iA and tending to infinity in the upper half-plane, avoiding the 
domain Di . Note that the union of contours C U Coo U C* U divides the complex plane in half. We say 
that the left (right) half-plane according to C U Coc U C* U is the half containing the negative (positive) 
real axis. For 3(A) > we then define 

L{X) + i6'"(A) , for A in the left half-plane according to C U Coo U C* U C^, 
L{X) - ie"{X) , for A in the right half-plane according to C U Coo U C* U C^. 
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(78) m := 



We note from H76|l that L{X) extends continuously to C and thus may be viewed as a function analytic for 
A € C+ \ Coo, where C+ denotes the upper half-plane. We introduce the notation 

(79) r(A) := 2e-^(^)/'^P(A) cos(6lO(A)/;i) . 

so the jump relation holding on A'o may be equivalently written in the form 

g-[g+(A)-g_(A)]/n q 

iT(A)e[2»QW+^W-9+W-f-Wl/'i el9+W-9-W]/ri 
The jump matrix in (|80|l can be factored as follows: 

ir(A)e[2^Q(^)+^(^)-9+(^)-9-(^)l/'i e[9+(A)-s-(A)]/?i 



(80) 



N+(A) = N_(A) 



XeKo. 



(81) 



T(A)"l/2 _^2^(-;^)-l/2g[-2iQ(A)-L(A)+2s_(A)]/R 

T(A)i/2 



jg[-2iQ(A)-L(A)+g+(A)+ff_(A)]/ft 
[2iQ(A)+"L(A)-g+(A)-s_(A)]/n q 



T(A)^/2 _j'7^(-;^)-l/2g[-2jQ(A)-L(A)+2g+(A)]/?i 

r(A)-i/2 

Such a factorization makes sense because we will see later (in § 13.7(1 that T(A) ^ 1 as iV ^ oo, so the 
fractional powers are well-defined for large enough N as having similar asymptotics. converging to 1 as 
N — > oo. The left-most (respectively right-most) matrix factor is the boundary value on Kq taken by a 
function analytic on the "minus" (respectively "plus") side of Kq. Let Kj^ denote a contour arc in D_i 
connecting the point Aq with the point A = —e, and oriented in the direction away from Aq. Let Dl denote 
the region enclosed by Kq, Kl, and the interval [— e,0]. We introduce a new unknown 0(A) based on the 
above factorization as follows: 

T{X)-^/^ _ir(A)-l/2e[-2iQ(A)-L(A)+2g(A)]/ri 

T(A)i/2 



(82) 



0(A) := N(A) 



(83) 



0(A) :=N(A) 



XeDL 



T(A)~l/2 iT(A)-l/2g[-2iQ(A)-L(A)+2g(A)]/?i 

r(A)i/2 

and elsewhere in the upper half-plane we set 0(A) := N(A). For A in the lower half-plane, we define 0(A) 
so as to preserve the symmetry 0(A*) = a20{X)*a2- 



Proposition 1. The matrix 0(A) has no jump discontinuity across the real intervals (— e,0) or (0, e), and 
thus may be viewed as an analytic function for A G C \ {K^i U Kq U Kl U Ki U /v* j U Kq U U K^) that 
takes continuous boundary values from each region where it is analytic. 



Proof. The boundary value taken by 0(A) on (0, e) from the upper half-plane is 

r(A)-i/2 



(84) 



0+(A)=N+(A) 



ir(A)"^/2g[-2jQ(A)-L(A)+»e«(A)+23(A)]/n 







9"(A)/S) 



T(X)-^/'^ . p(;^)g[-i(A)-i-»e°(A)]/n(g»e°(A)/;i g- 

which follows from H82|l . where we used the fact that (0, e) is in the right half-plane according to C U Coo U 
C* U to write L{X) in terms of L{X) with the help of (|78() . Here all of the quantities in the exponent 
are analytic functions on (0, e), and T{X)~^^'^ is interpreted in the sense of its boundary value taken on (0, e) 
from the upper half-plane. From the conjugation symmetry relations satisfied by 0(A) and N(A), it then 
follows that the boundary value taken by 0(A) on (0, e) from the lower half-plane is 

(T{X*)*)-^/^ ■ p(A)-ie[^(^)+*^°(^)l/'*(e'''°(^)/'* + e-»«°(^)/'^) 

-i(r(A*)*)-l/2g[2jQ(A)+L(A)+ieO(A)-2g(A)]/?i {T{X*)*)-'^^^ 
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(85) 0-(A)=N_(A) 



Here we have used the relations (|68|) . (|77|l . and H73|) to simpHfy the exponents. Since the matrix factor 
appearing in 185() has determinant one, to compute the jump relation for 0(A) across the interval (0, e) we 
will need to know the boundary value of the product T(A)~i/2(T(A*)*)~i/2 as A approaches (0, e) from the 
upper half-plane. First, note that for 3(A) > in the right half-plane according to C U Coo U C* U we 
may use H77|l and (|73|l to find that 

(86) T(A)r(A*)* = (1 + e2'»"(^)/'i)2 . 

(Similarly, if 3(A) > and A is in the left half-plane according to C U Coo U C* U C^, then the identity 

T(A)r(A*)* = (1 + e-2'e''(^)/'^)2 



(87) 

holds.) The pointwise asymptotic (sec 13. 7|) that r(A) 



-1/2 



1 as — !■ OO then gives that 



-1/2 ^ 



(1 



'{\)/hy 



(1 + e-2'e"(A)/R-)-i 



Ae (0,e), 
Ae (-6,0) 



Using this fact, one substitutes H85|) into (TSJ, and then H75II into H84|l . It is then an elementary calculation 
to deduce that 0+(A) ^ 0_(A) for A G (0,e). 

Starting with (|83|l , and proceeding in a similar way as we did to arrive at H84() , we find that the boundary 
value taken by 0(A) on the real interval (— e, 0) (which lies in the left half-plane according to CUCooUC*UC^) 
from the upper half-plane is 



(89) 0+(A) = N+(A) 
By conjugation symmetry, we then find 

(90) 0_(A)=N_(A) 







(r(A*)*)-i/2 p(A)-ie[^(^)-*''°(^)l/'*(e''''°(^)/'* + e-''''°(^)/'*) 

\)+L(A)-i6(''(A)-2g(A)]/S 

then shows that 0+(A) = 0_(A) for A £ (—£,0) as 



^(j^yyyl/2^[2tQ(X)+L(X)^te°(\)-2g{X)]/h (T(A*)*)~l/2 



Combining ^7^, and with the help of 

well. 

That the boundary values taken by 0(A) from each component of the complex plane where it is analytic 
are in fact continuous functions along the boundary even at self-intersection points follows from the fact that 
this was initially true for m(A), and is preserved by each of our substitutions to arrive at the matrix 0(A). 
(But it is also straightforward to verify this directly, even at the points A = — e, A = 0, and A = e.) □ 

The contours where 0(A) has discontinuities in the complex plane are illustrated with black curves in 
Figure 

3.4. The choice of g{X). Bands and gaps. Up until this point, the contours have been more or less 
arbitrary, with the only conditions being that the region D^i contain the imaginary interval [0,i^] and 
the branch cut C. Likewise, the function g{X) remains undetermined aside from its analyticity properties 
relative to the contours and the symmetry property (|^ . We now will describe how to choose both the 
system of contours and the function g{X) to render the construction of 0(A) asymptotically tractable in the 
limit N —f oo. 

Recall that on the contour Kq we have the jump relation H8U|I for N(A). Since Kq lies in the right-half 
plane according to C U Coo U C* U C^, wc may write the corresponding jump relation for 0(A) in the form 
(91) 

Q jg-[2iQ(A)-|-L(A)-je°(A)-s+(A)-g_(A)]/?i 
^p[2tQ{X)+L{X)-te°{X)-g+iX)-g^{X)]/h q 



0+(A) = 0_(A) 



XeKn 



Similarly, the jump relation satisfied by 0(A) for A on the contour A'_i follows from (|70|) and can be written 
in the form 

e-[s+(A)-9-(A)]/rx 
iS'(A)e[2»Q(^)+-f^(^)-»«°(^)-f+(-^)-9-(^)l/'i elg+W-g-W]/ri 
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(92) 



0+(A)-0_(A) 



A e , 




Figure 6. The contours Ki, Kq, K-i, and K^, and their conjugates are the contours of 
discontinuity of 0{X) and are shown in black. Superimposed in blue are the curves C and 



where we define 

(93) Six) := P{\)e-^^^^/''\ 

The key observation at this point is that whether A G Kq or A G A'_i, the exponents of the entries in the 
jump matrix have the same form. 

We suppose that the contour loop Kq U K-i across which g{X) is allowed discontinuities in the upper 
half-plane is divided into two complementary systems of sub-arcs called bands and gaps. The characteristics 
of bands and gaps are most easily phrased in terms of two auxiliary functions defined along Kq U K^i that 
are related to the boundary values of 5(A): 

(94) 0(A) :=*(g+(A)-.9_(A)), 

(95) <Pi\) := 2iQ{\) + L(A) - iO^iX) - g+{X) - g^{X) . 
In terms of these functions, we have the following definitions. 
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• A band is an arc along which the following two conditions hold: 

(96) 't'W = imaginary constant, 0(A) = real decreasing function. 

We will require that the arc Kq is a band, along (possibly) with some sub-arcs of K-i. 

• A gap is an arc along which the following two conditions hold: 

(97) 3?(0(A)) < , 61(A) = real constant. 

We will require that the terminal arc of A'_i (that meets the real axis at the point A = — e) is a gap, 
along (possibly) with some other sub- arcs of K-i. 
By "decreasing" in (|96|l we mean in the direction of orientation, namely counterclockwise. Note however, 
that while the orientation of a contour containing a band is essentially an arbitrary choice, the condition 
that 9{X) be decreasing has intrinsic meaning because 9{X) is by definition proportional to .9+ (A) — g-(A) 
which changes sign upon reversal of orientation. 

Given an arbitrary set of contours, it may be the case that for no function .g(A) defined relative to 
these contours can the conditions H96|) and H97|l both be satisfied by choice of systems of bands and gaps. 
Specification of the contours is thus part of the problem of finding g{X). 

To find g{X) and the contours along which its discontinuities occur, we suppose that the number of bands 
along Kq U A'_i is known in advance, and show that certain conditions are then necessary for the existence 
of g{X) for which the conditions H9()|) and H97|l both hold. Taking a derivative, we see that the function g'{X) 
must satisfy the following four conditions 

(98) g'+{X) + g'_{X) = 2iQ'{X) + L'{X) - ie^'{X) , for A in the bands, 

(99) 5+ (A) - g'_(A) = , for A in the gaps, 

(100) 5'(A)+5'(A*)* = 0, 
and 

(101) ,9'(A) = 0(A-2), as A ^00. 

(Here we are assuming that differentiation with respect to A commutes with the limit of taking boundary 
values. That this is justified will be seen shortly.) These four conditions amount to a scalar Riemann-Hilbert 
problem for the function g'{X). 

To solve for 5' (A), it is most convenient to first modify the unknown to eliminate the L'{X) term from 
(jnU. Let Cf denote the contour Kq U K-i U (— e, 0) U C oriented in the direction from A = to A = iA. Set 

(102) /c/(A):=^ / XeC\Cf, 
and then 

(103) J{X):= fu{X)- !u{X*)\ XeC\{CfUC}). 
The change of variables we make is 

(104) h{X) := g'{X) + /(A) . 

Note that if A ^ then since Kq U A'_i U (— e, 0) = dD^i, we may write /c/(A) in the equivalent form 

(105) MA) = i/A = irA. A^i..,, 



Jc rj- X 2 Jo 7/ - A 

from which it is clear that /[/(A) = A/2X + 0{X^^) as X^ 00. Therefore /(A) = 0(A~^) in this hmit. Also, 
from (|103|l . we have 

(106) /(A) + /(A*)* =0. 
The formula H105|l also shows that 

(107) f{X)^-l-L'{X), XeC\{D^,UD*_,), 



2 
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and from the Plcmclj formula, at each point of Kq U K^i and C we have 

(108) /+(A) - /_(A) = -11 = i0"'{X) , A e Ko U K-i U C , 

while from the condition H1U6|) and (|77|) we get that /(A) has a jump discontinuity across the real interval 
(— e, 0) given by 

(109) /+(A)-/_(A) = -27r, Ae(-e,0), 

where orientation from left to right is understood (so /+ is a boundary value from the upper half-plane). It 
follows that h{X) is a function analytic for A G C \ A'o U U C U Kq U Kl^ U C* U (— e, 0) whose boundary 
values are related to those of g'{X) as follows: 

(110) h+{X)+h^{\)=g'_^{\)+g'_i\)-L'{X)+i9°'{X), XeKo^K-i, 

(111) /i+(A)-/i_(A)=,g;(A)-,g^(A)-7r, A e i^To U A'_i U C , 
and 

(112) h+{X) - h^{X) = -2t: , Ae(-e,0). 
Consequently, h{X) is characterized by (|112|l and the following four conditions: 

(113) h+{X) + h_{X) = 2iQ'{X) , for A in the bands of Kq U A'_i, 

(114) h+{X) — h-{X) = — TT , for A in the gaps of Kq U A'_i and on C, 

(115) h{X) + h{X*)* = , 
and 

(116) h{X) = 0{X-^), asX~^oo. 

These conditions make up a scalar Riemann-Hilbert problem for h(X). Note that the boundary values taken 
by h{X) are continuous except at A = 0, where a logarithmic singularity is allowed (and necessary, to cancel 
the corresponding singularity of /(A)). 

To solve for h{X) and thus obtain g'{X), wc suppose that the cndpoints of the bands along Kq U A'_i are 
given: 

(117) Ao, Ai, . . . , Ag , in order along Kq U Ar_i, for G even. 

(Despite similar notation, these are not directly related to the soliton eigenvalues Ajv.fc.) We then introduce 
the square-root function i?(A) defined by the equation 

G 

(118) R{xr = Y[{x-x,,){x-x:), 

and the condition that the branch cuts are the bands of A'o U K^i and their complex conjugates, and that 
R{X) = X'^^^ + 0{X^) as A ^ oo. We attempt to solve for h{X) by writing it in the form 

(119) h{X) = R{X)k{X) , 

for some other unknown function fc(A). Since R{X) = R{X*)* and since the jump discontinuities of R{X) 
are restricted to the bands and their conjugates, where they satisfy R+{X) = — i?_(A) for A in the bands of 
Kq U Ar_i, the conditions imposed on h{X) take the form of conditions on fc(A) as follows: 

(120) fc+(A) - fc_ (A) = ^^-YT , for A in the bands of A'o U K_i, 

TT 

(121) fc+(A) — fc_(A) = — , for A in the gaps of A'o U A'_i and on C, 

R[X) 
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(122) 



k+{\)-k^{X) = -^^, for A e (-6,0), 



(123) k{X) + k{X*)* = , 
and 

(124) A:(A) =0(A~('^+3'), as A -> cx). 

We allow k{X) to have singularities at the band endpoints. as long as the product R{X)k{X) is regular there, 
and to have a logarithmic singularity at A = 0. 

Since the differences of boundary values of fc(A) are known, it is easy to see that it is necessary given this 
information that k{X) has the following form: 
(125) 



ku{X):=klj'>{X)+k'^'>{X), kl^'iX):^- " ' , , , ^^(A) , ^, . , , 

and 

(126) k{X):=kuiX)~ku{X*)*. 
Because a residue calculation gives 

(127) fc«(A) - klj\xr = ^ - 2^t5a,o , 
the formula for fc(A) becomes 

^ _0,VA„.^ 1.(2) 1.(2) r 

R{X) 

The required properties (|12U|l - p23|) are satisfied by this expression, and its singularities are easily seen to 
be of the required types. Here we can see also that as a consequence of the analyticity of the densities of the 
Cauchy integrals used to define fc(A) and hence g{X), differentiation of ^(A) clearly commutes with taking 
boundary values, at least away from the endpoints of the bands. 

On the other hand, the decay condition (|124|l is not necessarily satisfied by our formula for fc(A). By 
explicit expansion of (|128|l for large A, we see that for (I124II to hold the following additional conditions are 
required. For any integer p > 0, define 



(128) k{X) = - 2it5Gfi + k^J' (A) - k'^J> (A*)* 



(129) Mp(Ao, 



rf drj 

C/\bands 

R{ri) 



If G = 0, then the condition (|124|l requires that 

A/o(Ao) = x + 2tao, 

(130) . ^ 

Afi(Ao) = xao + 2t[al- -bl 



while for even integers G > 2, H124|l requires that 

Mp(Ao,...,Ag) = 0, forO<p<G-2, 

A'/g-i(Ao,...,Ag) - 2t, 

G 



(131) AfG(Ao,...,AG) = :E + 2t^a„, 

G G / ^ \ G G 

Mg+i(Ao,...,Ag) = a„ + 2i ^ f - -6^ j + 2t ^ ^ a^a^ 



Ti=0 n— ^ n— m— n+l 
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In these formulae, a„ := 5R(A„) and 6„ := 3(A„). These are G + 2 real constramts that must be satisfied by 
choice of the G + 1 complex numbers Aq, . . . , Xg in the upper half-plane. 

For each configuration of contours and band cndpoints consistent with the conditions ()13U|I or we 
therefore obtain a candidate for the function g{X) by the formula 

POO POO 

(132) 5(A) - / g'(r;) dr? = - / [R{v)k{v) - /(r/)] dry , 

where the integration path is an arbitrary path from A to infinity in the region C \ (£*-! U -D* i) if A lies 
in this region as well, whereas if A e -D-i U D*_i, then a first component of the path lies in this region and 
connects A to — e, followed by a second component that coincides with the real half-line (— oo, — e). We may 
then attempt to enforce on g{\) the conditions that 3fJ((/)(A)) — within each band and 3(6'(A)) = within 
each gap. For ?i = 1, . . . , G/2, let A„ denote a simple closed contour with positive orientation that surrounds 
the band with cndpoints \2n-1 and A2,i and no other discontinuities of g'(A). The definition H94|l and the 
integral formula (|132|l for g{\) shows that 6{\) = in the terminal gap of K_i (from the point \ = \q 
along A'_i to the point A = — e). Therefore, 9{\), already made constant in the remaining gaps via the jump 
conditions imposed on .g'(A), will have a purely real value in each gap if 



(133) 



0, 



n = l,...,G/2. 



Similarly, for n ~ I, . . . , G/2, let r„ be a contour arc representing the gap in A'_i between A2T1-2 and A2,i-i. 
The definition H95|l shows that (due to the symmetries H77|) and (I68II ') 0(A) is already purely imaginary in 
the band Kq. Consequently, its constant value in the remaining bands will be purely imaginary also if 



(134) 



[2g'(77) - 2iQ'{ij) - L'irj) + tO'^'iv)] diA ^ 



n = 1, . 



,G/2. 



If G is an even positive number, then the conditions (|131|1 . (|133|1 . and H134|l taken together are real equations 
sufficient in number to determine the real and imaginary parts of the complex cndpoints Ao, . . . , Xg along 
Kq U K-i. If G = 0, then there are no conditions of the form (|133|1 or H134|l and the equations (|130|l are 
expected to determine the single endpoint Aq. 

The conditions (|133|l and H134|) can be written in a common form that is useful for computation. First, 
note that since 5+(?7) — g'-ii]) = g'iv) for V any gap r„, from H110() wc get 



(135) 



2.g'(r,) - 2zQ'(7y) - L'(r;) + i9"' {rj) = h+{rj) + /i_(r?) - 2tQ'{Tj) 



2iQ'{rj) 



Riv) 

■fcf (77) + fci2)(r?)l , 



?7 e r„ , 



Then, using (|119|l . 

(136) 2g'(77) - 2iQ'{rj) - L'{i^) + lO"' {r^) = R{rj) 

Using H128|l in the case G > 0, we then get 

(137) 2g'{Tj) - 2iQ'{ij) - L'{ri) + (77) = i?(r?) 

where fc(2)(,y) -.^ fc(2) ^ \^ _ 1.(2) 
shows that k^^\rf) + k^_^' {vi) = ^^(77), where 

(138) Y{i^):^Yu{ii)^Yu{v*)\ and Yuiv) ■ , ■ w w , 

^ J -00 R{v){v - v) JcR{^){^-v) 

(Note that the integrand in Yij(ri) has a jump discontinuity on the real axis at = due to the factor R{ri) 
in the denominator.) To simplify the conditions 1133|l . we let /„ denote the band arc connecting A2„_i to 
A2„, and note that H133|l can be written in the equivalent form 



jj (A) — fcjy (A*)*. For 77 in a gap r„ of 1 , an elementary contour deformation argument 



dv 



(139) 



1, 



,G/2. 



Indeed, this form is more natural given the connection of these conditions to the net change in the function 
6{X) as A moves through the band /„. Using (|lll(l . wc have 

(140) g'_^i,j)-g'_iv) = h+{n)-h^i7^)+n, 77 € /„ , 
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and then from H119|) . 



(141) 



g'+iv) - g'-iv) = R+iv) 



k+iv) + k-{r/) 



R+iv). 

because R(r]) changes sign across the branch cut /„. Another contour deformation argument then shows 
that the combination fc+(?y) + fc-(fy) + T^/R+iv) ms-Yi for V S Im be identified with the same function Y{j]) 
as defined by (|138|l . Therefore, the conditions H133() and 1] 134(1 may be written similarly as 



(142) 



i?„(Ao, . . . , Ag) 



K(Ao,...,Ag) 




R{-q)Y{'q) drj 



= 



l,...,G/2, 



where the paths of integration lie in the region of analyticity of the integrand. Now, strictly speaking, this 
does not amount to a definition of functions i?„(Ao, . . . , Ag) and 14,(Ao, . . . , Ag) because the integrals are not 
individually independent of path due to monodromy about the branch cuts of R{ri). However, the totality 
of the conditions (|142|l is clearly independent of any particular choice of paths (for example, adding to the 
path from Aq to Ai a circuit about the branch cut of Rirj) connecting Ai and A2 amounts to adding to Vi a 
multiple of which is zero on a configuration satisfying H142|l ). 

3.4.1. Whitham equations. The endpoints Ao,...,Ag are determined implicitly as functions of x and t 
through the equations En{v) ~ for 71 = 1, ... , 2G+2, where the unknowns are v = (Aq, . . . , Ag, Ag, . . . , A^)"^ 
and the equations are 

(143) 

and 

(144) 

for n in the range n = 1, . 
(145) 



En{v) := V(„+i)/2(u) , for n odd, 

En{v) Rn/2{v) , for n even, 
. , G, and then 

E„{v) := M„„G-i(v) , forn = G + 1, . . . , 2G - 1, 



(146) 



E2Giv) := Mg-i{v) - 2t , 



(147) 



and 



(148) 



E' 



■2G+2 



E- 



2G+1 



{v) Mg{v) 



2t y ^ a„ 



n=0 



(v) ■.= Mg+i{v) 



■ G 



2t 



G 

E 

ri=0 



2t 



G G 

E E 

n— m— n+1 



In these equations, = {Xk + A^.)/2, = (Xk — AJ!)/(2i), and A/„(t/), i?„(w), and Vn{v) stand for the 
complexification of the corresponding real quantities. That is, 

1 f 7f d7] 1 r 7f d7] 1 rf dri iWP+^ r e*(P+i)^ dO iWP+^ f° e*(P+i)^ d9 

2 Jclm'^ Jc'RM^^ J-wW)^~ 



(149) Mpiv) .- ^ , ^^^^ ^ , ^^^^ . ^ , ^^^^^ , ^ 1^ R{WeJ^) 4 J_, RiWe^^) 

where the contour G* is taken to be oriented from rj = —iA to = and Vt^ > is a sufficiently large 
number (if p < G then we may pass to the limit W ^ 00 and drop the last two integrals). The complexified 
Mp agrees with the expression defined by 1)129(1 when A^ = + ibk and AJ! = — ibk with and bk being 
restricted to real values. The complexified Mp is a function of the independent complex variables Xk and A^ 
through the branch points of the function R in the integrand. Similarly, 

(150) Rn{v) i / " Riv)Y{v) dv - I f R{v)Y{7i) dv , 
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and 



(151) 



-^2ti — 2 



R{v)Y{ii)dji-- 



where in both cases the integrals in the two terms are taken over complex-conjugated paths. Simple contour 
deformations near the endpoints of integration show that in each case derivatives of and Vn{v) with 

respect to any of the Vk at all can be calculated by differentiation under the integral sign (even if Vk is one of 
the endpoints of integration). Once again, these complexified quantities, while functions of the independent 
complex variables vi,. . . ,W2G+2, agree with the previous definitions when Vk+c+i = v^. Of course we are 
only interested in those solutions of the equations En{v) = that have this conjugation symmetry. 

If v{x,t) is a differentiable solution of the equations En{v) = for n = 1, . . . , 2G+2, then we may calculate 
dv/dx and dv/dt by implicit differentiation. Thus, 



(152) 



dE dv _ BE ^^^^ dE dv _ dE 
dv dx dx ' dv dt dt ' 

where dE/dv denotes the Jacobian matrix of the E„ with respect to the Vk holding x and t fixed, while 
dE/dx and dE/dt are the corresponding vectors of partial derivatives with respect to x and t (holding the 
Vk fixed) . Clearly, these latter partial derivatives contain no explicit x and t dependence because the En are 
all linear functions of x and t. In fact, it turns out that when v = v{x,t) satisfies the equations E{v) = 0, 
the Jacobian matrix can also be expressed up to a diagonal factor in terms of v alone (x and t may be 
eliminated). Indeed, direct calculations show that 



(153) 



and 



(154) 



ovk 2 ?7 - Vk 



dMo 
dvk 



(155) 



dE„ 



dvk dvk 
Also by direct calculation 

(156) 

so on a solution v{x, t), 



^" RiTj)d7J 



A! 



V - Vk 



= ■ Jnk{v) 



dM< 



dvk 



for 1< 71 < G . 



dEn _ 1 ^ 



Vk- 



dE„- 



dvk 



for G + 2 < n < 2G + 2 . 



(157) 



9 En _ n^idEa+i 



V,. — — =: Jnk(v)- 



for G + 1 < n < 2G + 2 , 



dvk ovk " ovk ovk 

assuming that G > 0. This proves that on a solution v = v{x,t), the Jacobian matrix can be expressed in 
the form 



(158) 



dE (dMo 
-=J(.)-diagf^, 



dM, 



' dv2G+2 , 

where the matrix 3{v) is an explicit function of v whose elements arc defined by ()155ll and H157() . Its 
determinant is nonzero as long as the Vk are distinct. Applying Cramer's rule to 1)152(1 . wc then find that 



for fc = 1, , 
det J*''* 



(159) _£+Cfc(t;)^ = 0, for fc = l,...,2G + 2, 
where 

(160) cki^)- ^^^j,,., 

where J'^'* (respectively J'^'^) denotes the matrix J with the fcth column replaced by dE/dt (respectively 
dE/dx). The system of quasilinear partial differential equations (|159|l satisfied by the endpoints v{x,t) 
is automatically in Ricmann invariant (diagonal) form, regardless of how large G is. These equations are 
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frequently called the Whitham equations. They clearly play a secondary role in our analysis, as they were 
derived from the algebraic equations En (v) = which are more fundamental (and in particular encode initial 
data). 

3.4.2. Inequalities and topological conditions. Subject to being able to solve for the cndpoints Ao,...,Ag, 
there is now a candidate for the function g{X) associated with each even nonnegative integer G. We refer to 
such a guess for g(X) below as a genus-G ansatz. The selection principle for the number G is that it must 
be chosen so that the inequalities 

< for A in each gap of Ki^ U K^i, 

(161) 

6{\) is strictly decreasing along bands of Kq U K_i 

both hold true. Given that the actual band and gap arcs have not yet been chosen, these conditions are 
really topological conditions on the level curves of the real part of the integral 

(162) /(A) = / i?(r;)r(77)dr;. 



(These level curves are also known in the literature as the orthogonal trajectories of the quadratic differential 
i?(A)^y(A)^(iA^.) A band arc of i^o U K^i must coincide with a level curve of 5R(/(A)) connecting the origin 
with Ao, or \2n-1 with X2n for n ~ 1, . . . , G/2. Furthermore, it must be possible to choose the remaining 
arcs (gaps) so that they lie in the region where 5R(/(A)) is less than at either endpoint. 



(163) 



0+(A) = 0_(A) 



3.5. Third modification: opening lenses around bands (steepest descent). In terms of the functions 
9{X) and 0(A), the jump discontinuity of 0(A) across A'„i takes the form 

^5(A)e*(^)/'^ e-'^i^y^ 

Assuming that K-i remains bounded away from the imaginary interval [0,i^] of accumulation of poles for 
m(A) by a fixed distance, we have by a midpoint rule analysis for Riemann sums that 5'(A) = 1 + 0{K). If 
A is a point in a gap Tj C K_i, then 9{\) = 6j € M, and 3?((/)(A)) < 0, so the jump matrix in H163|) is an 
exponentially small perturbation of the constant (with respect to A € Tj) jump matrix e**^"'^/''. 

On the other hand, in direct analogy with the factorization l|81f) . the jump discontinuity of 0(A) across 
K-i can be written in factorized form: 



iS{X)e 



4>{\)/h 







(164) 



5(A)-i/2 




-i5'(A)"i/2g-0(A)/?igje(A)/ri 
5(A)i/2 



ie-'t'WIh 



5(A)i/2 




-i5(A)-i/2e-0(A)/?.g-*e(A)/rx 



5(A)- 



1/2 



This factorization is useful for A in a band Ij C i^-i. (Recall that we are also assuming that the contour 
Kq is itself a band in its entirety; thus Kq = Iq). Indeed, let £1^ (respectively ilj) be a lens-shaped domain 
lying to the left (respectively right) of the band Ij C Let iKj be the purely imaginary constant value 

of 0(A) in the band Ij. We introduce a new unknown P(A) based on this factorization as follows: 



(165) 



P(A) := 0(A) 



5(A)-i/2 




iS{X) 



^iji -iKjlh -iB{X)lh 



5(A) 



1/2 



Aef)+, 



(166) 



P(A) := 0(A) 



5(A)-i/2 




-i5(A)"^/2e"'''j/'^e'*(^'/'^ 



5(A) 



1/2 
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for all other A in the upper half-plane where 0(A) takes a definite value we set P(A) = 0(A), and finally 
for all A in the lower half-plane we set P(A) = (T2P(A*)*(T2- In writing down this change of variables we are 
making use of the fact, apparent from the explicit formula for .g'(A). that the function 6{X) has an analytic 
continuation from each band Ij C A'_i to the regions fi^. The contours in the complex plane across which 
P(A) has jump discontinuities arc shown with black curves in Figurc|3 A simple Cauchy-Riemann argument 




Figure 7. The discontinuity contours for P (X) . 

taking into account the monotonicity of the real-valued analytic function 0{X) in the bands then shows that 
on the boundary contours dfl^ (not including the band Ij, as shown in Figure0l the jump matrix converges 
to the identity matrix as ?i — > 0. The convergence is uniform away from the endpoints of the band, with a 
rate of convergence 0{h). 

These heuristic arguments will be used to suggest in 13.61 a model for P(A) we call a parametrix, and 
then they will be recycled in 13.71 to prove that the parametrix is indeed an accurate model for P(A). For 
now, it suffices to note that the matrix P(A) is the unique solution of a matrix Ricmann-Hilbcrt problem 
that by our explicit steps is equivalent to the discrete Ricmann-Hilbcrt problem satisfied by m(A). This 
problem is the following. Seek a 2 x 2 matrix P(A) with entries that are piecewise analytic functions of A in 
the complement of the contours if-i, Kq, Ki, Kl, the lens boundaries d^f, and their complex conjugates 
such that 
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(167) 



(168) 



(169) 



(170) 



(171) 



The boundary values taken on each arc of the discontinuity contour are continuous along the arc and 
have continuous extensions to the arc endpoints. The boundary values are related by the following 
jump conditions. For X £ Ij C K-i or for A S /q = -fCo, 



P+(A)=P_(A) 
On the arc Kl we have 
P+(A)=P_(A) 





.IKj j h 



iKj / h 



A e /, 



r(A)-i/2 j2.(;^)-i/2g-,«o/?^e-»e(A)/rx 
r(A)i/2 



where 9{X) refers to the value of the function 6{X) analytically continued from Iq to K^. On the arc 



P+(A) = P-(A) 



T(A)i/2 iT(A)~^/2e~*""/'*e*''(^)/'^ 
T(A)-i/2 



1 

iS'(A)e*«o/'*e'(2e''(^)-e(^))/'i i 



A € Ki 



where again 0(X) refers to the value analytically continued from Iq. On a gap Tj C K-i, 



P+(A)=P_(A) 



iS'(A)e'^(^)/'^ e^'"^/'' 



X e r, 



(This holds with 9j = on the final gap of K^i from Ag to — e because g{X) is analytic in this gap 
so 9{X) = 0.) On the lens boundaries dft^, 

5(A)Tl/2 ,5(^)-l/2g-,:.,/;igT*e(A) 

5(A)±i/2 



P+(A) = P_(A) 



A e dilf , 



(172) 



where here 9{X) refers to the value analytically continued from . Finally, the jump relations satisfied 
in the lower half-plane are by definition consistent with the symmetry P(A) = (T2P(A*)*cr2. 
The matrix P(A) is normalized so that 

lim P(A) = I. 

A— ^oo 



The function defined in terms of P(A) by the hmit 



(173) 



ipN{x,t) ■= 2i lim APi2(A) 

A — >oo 



is the A^-soliton solution of the semiclassically scaled focusing nonlinear Schrodinger equation. 

3.6. Parametrix construction. Building a model for P(A) consists of two steps: (i) dealing first with the 
asymptotic behavior of the jump matrix away from the endpoints of the bands, and (ii) local analysis near 
the endpoints. Both of these constructions have been described in detail elsewhere, and we will give just an 
outline of the calculations. 



3.6.1. Pointwise asymptotics: the outer model problem. As iV ^ oo (so in particular h ~ Hn 0) the jump 
matrix defining the ratio of boundary values in the Riemann-Hilbcrt problem for P(A) may be approximated 
by a piecewise constant (with respect to A) jump matrix that differs from the identity matrix only in the 
bands of Kq U K-i and their complex conjugates and also in the nonterminal gaps of K-i. Thus we may 
pose another Riemann-Hilbert problem whose solution we hope we can prove is a good approximation in a 
certain sense of P(A). We seek a 2 x 2 matrix function P(A) that is piecewise analytic in the complement of 
the bands and nonterminal gaps and their complex conjugates such that 

• The boundary values taken on each band or gap of the discontinuity contour arc continuous along 
the arc and have singularities of at worst inverse fourth- root type at the band/gap endpoints. The 
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boundary values are related by the following jump conditions. For X E Ij (Z K-i or for A e /o = A'o, 



(174) P+(A)=P_(A) 
For A in a nonterminal gap Fj C K^i 

(175) P+(A)=P_(A) 



ie-^'^o/f' 



g-ifj/S 



A G L 



A e F,- 



Finally, the jump relations satisfied by P(A) in the lower half- plane are by definition consistent with 
the symmetry P(A) = cr2P(A*)*cr2. 
• The matrix P(A) is normalized so that 

(176) lim P(A) =1. 

A — voo 

This Riemann-Hilbert problem can be solved by first introducing an auxiliary scalar Riemann-Hilbert prob- 
lem with the aim of removing the jump discontinuities of P(A) across the nonterminal gaps as expressed 
by the jump conditions (|175|) while converting the jump matrix for P(A) in all of the bands to the matrix 
CTi . The two columns of the resulting matrix unknown thus may be considered to be restrictions of a single 
vector-valued analytic function on a hyperelliptic Riemann surface X constructed by identifying two copies 
of the complex A-plane across a system of cuts made in the bands and their complex conjugates. This Rie- 
mann surface has genus G, which explains our terminology for the "genus" of a configuration of endpoints 
for the g-function. The vector-valued function defined on the surface X that leads to the solution of the 
Riemann-Hilbert problem for P(A) is known as a Baker-Akhiezer function. It can be expressed explicitly in 
terms of the Riemann theta function of X, as can P(A). The details of this construction can be found in 

nni. 

Two features of this solution are important for the subsequent steps in the analysis. First, the dependence 
of the solution on h enters through the real quantities 9j/h and Kj/h, which determine a point in the real 
part of the Jacobian variety of X, topologically a torus of dimension G. Essentially, this point is a phase 
shift in the argument of the Riemann theta functions used to construct the solution, and in particular this 
implies that the matrix P(A) remains uniformly bounded for A away from the band/gap endpoints in the 
limit ft — > 0, even though the phase point oscillates wildly in the Jacobian in this limit. Next, for A near the 
band/gap endpoints, the matrix P(A) exhibits a singularity of a universal type that, while a poor model of 
P(A) near the endpoints, nonetheless turns out to match well onto another matrix function that is a better 
model. 



3.6.2. Endpoint asymptotics: the Airy function local parametrix. To determine what this better model should 
be, it suffices to fix a sufficiently small neighborhood Uk of each band/gap endpoint A^, and to find a matrix 
that exactly satisfies the jump conditions of P(A) in this neighborhood. Such a matrix can be found because 
the jump matrices restricted to Uk can be written in a canonical form with the use of an appropriate conformal 
mapping (Langer transformation) taking Uk to a neighborhood of the origin. Once the jump matrices in 
Uk are exhibited in canonical form, a piecewise analytic matrix function satisfying the corresponding jump 
conditions can be written down explicitly in terms of Airy functions. Next one observes that in fact there 
are many piecewise analytic matrices defined in Uk satisfying the exact jump conditions of P(A), all differing 
only by multiplication on the left by a matrix factor analytic in Uk- The choice of this factor can be used 
to single out a particular local solution that is a good match onto the explicit matrix P(A) on the boundary 
dUk of Uk- Specifically, one chooses the factor so that the resulting local solution, which we will call Pfe(A), 
satisfies 

(177) Pfc(A)P(A)-i =I + 0(ft), 

as fi, ^ 0, uniformly for A € dUk- 

The Airy function local parametrix is described in detail, for example, in PP. Here, we need a slight 
modification of the construction of jj, because the jump matrices for P(A) restricted to Uk involve the 
function S'(A), and also the function r(A) in the case of t/o- We can easily remove these functions from the 
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(178) 



Q(A) = P(A) 



jump matrices by making a local change of variables in Uk as the first step in the construction of Pfc(A). 
Suppose first that fc > and in the part of Uk lying outside of the lenses we set 

5(A)-i/2 
5(A)i/2 

while in the rest of Uk we set Q(A) = P(A). It is then easy to check that the matrix Q(A) satisfies the same 
jump conditions as does P(A) but with 5(A) simply replaced by 1. This turns out to be a near-identity 
transformation since S{X) = 1 + 0{h) uniformly for X Cz Uk- Next, consider the jump conditions satisfied by 
P(A) in Uq. In the part of Uq common to _D_i but outside of the lenses Di and Dl we set 

r(A)-i/2 
r(A)i/2 

while in the part of Uq outside both the lenses and D_i wc set 



(179) 



Q(A) = P(A) 



(180) 



Q(A) = P(A) 





T(A)i/2 



1 Oil" T(A)-i/2 

-i5(A)e''"i/'*e'(2e°(^)-e(^))/'i i 

and in the remaining parts of Uq we set Q(A) = P(A). Using the relationship between L{X) and L{X) valid 
in the right half-plane according to CU Coo U C* U it then follows that on TiDUq, KlDUq, and /q n C/q, 
the jump conditions satisfied by Q(A) are of the same form as those satisfied by P(A) but with S'(A) and 
T(A) both replaced by 1, while for Xe KiD Uq, 



(181) 



Q+(A) = Q_(A) 



From this point, the construction follows that in precisely, with Q(A) being studied within each Uk by 
means of an appropriate Langer transformation. 

3.6.3. Global parametrix. We now propose the following global parametrix, P(A), as a model for P(A) uni- 
formly valid in the whole complex plane. The matrix is well-defined globally with the exception of certain 
contours on which continuous boundary values are taken from each side: 



(182) 



P(A):=Pfc(A), for Ae C/fc, fc = G,...,G, 



(183) P(A) :=CT2P(A*)V2, for A e C/^ = 0, . . . , G, 
and 

(184) P(A) := P(A) , for A outside all neighborhoods Uk and their conjugates. 
3.7. Error analysis. We now argue that E(A) := P(A)P(A)^^ satisfies 



(185) 



lim A(E(A) - I) = 0{h) 



The basic properties of the matrix E(A) follow on the one hand from the conditions of the Riemann-Hilbert 
problem satisfied by the factor P(A) and on the other from our explicit knowledge of the global parametrix 
P(A). Clearly, E(A) is a piecewise analytic matrix in the complex A-plane. satisfying 



(186) 



lim E(A) =1, 

A — 'oo 



with jump discontinuities across the following contours: 

• Across the boundaries dUk of neighborhoods of endpoints in the upper half-plane, taken with coun- 
terclockwise orientation, we have 

(187) E+(A) = E_(A)P(A)P;,(A)-i, X e dUk . 

• Across the bands and nonterminal gaps of Ko U K-i outside the neighborhoods Uk in the upper 
half-plane, where both P and P = P have jump discontinuities, 

(188) E+(A) = E_(A)P_(A)v(A)v(A)-ip_(A)-i , 

where across the same contour P+(A) = P_(A)v(A) and P+(A) = P_(A)v(A). 
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• Across the portions of Ki, Kl, the terminal gap of K-i, and the lens boundaries 917 j. that he outside 
of the neighborhoods Uk in the upper half-plane, where only the factor P(A) is discontinuous, we 
have 

(189) E+(A) =E_(A)P(A)v(A)P(A)^S 

where the jump matrix v(A) is defined by P+(A) = P_(A)v(A). 

The jump discontinuities of E(A) in the lower half-plane are consistent with the symmetry E(A) = (T2E(A*)*cr2- 
In particular, E(A) is an analytic function inside all neighborhoods Uk of endpoints and their complex con- 
jugates because Pfc(A) is chosen to satisfy the jump conditions of P(A) exactly within Uk- 

This information means that E(A) itself is the solution of a matrix Riemann-Hilbert problem with given 
data. By exploiting a well-known connection with systems of singular integral equations with Cauchy-type 
kernels it suffices to estimate the uniform difference between the ratio of boundary values E_(A)~^E+(A) 
and the identity matrix I. In fact, wc will show that E_(A)^^E+(A) — I = 0{h) holds uniformly on the 
?i-indcpendcnt contour of discontinuity for E(A). From this estimate, the estimate (I185II follows from the 
connection to integral equations. 

To show that E_(A)^^E+(A) = I + 0{h) requires only a little more than the properties of g(A) and the 
global parametrix P(A) already established. We also need the asymptotic behavior of the functions S'(A) 
and T{X). Analogous functions are analyzed carefully in PP, so we just quote the results: 

• The function S{X) is analytic for A in the upper half-plane outside the region bounded by the curve 
C and the imaginary interval with the same endpoints [0, iA]. Uniformly on compact subsets of the 
open domain of analyticity we have S{X) = 1 + 0{h). 

• The function r(A) is analytic for A in the upper half-plane outside the region bounded by the 
curve Coo and the imaginary axis above lA. Uniformly on compact subsets of the open domain of 
analyticity we have T(A) = 1 + 0{h). 

These facts are enough to prove that E_(A)~^E+(A) = I + 0{h) on all contours with the exception of Ki 
restricted to a neighborhood of A = e and restricted to a neighborhood of A = — e. The jump matrix for 
P(A) on Kl involves both e^^W/f' and also e^^'^ The former is exponentially small as ft | by a 

Cauchy-Riemann argument for Ki sufficiently close to Kq] the choice of a sufficiently small but positive e is 
crucial to provide the decay where Ki meets the real axis. The latter is also exponentially small on the parts 
of Kl that are bounded away from the imaginary axis, because ^{2i0^{\)) = — 27r3fi(A), which dominates 
^^(^^(A)) for Kl close to Kq. But it is not immediately clear that an e > can be found so that the inequality 
3?(2i0°(A) — i9{\)) < persists along Ki to the real axis. However, taking a hmit of g+(A) — g'_{X) as A — > 
along Kq shows that 

(190) n{~ie'{0))=TT, 

so 3ff(2i6'°(A) — i0{X)) = — 7r5R(A) + 0(|Ap), which means that the inequality persists along Ki to A = e > 0, 
for e sufficiently small. A similar explicit calculation involving g'{X) near A = shows also that 3ff((^(A)) 
is decreasing linearly away from the origin along the negative real axis, which proves that while the limit 
of 0(A) as A approaches the origin along K-i is purely imaginary, the inequality 5R((/)(A)) < is satisfied 
strictly throughout the terminal gap of K^i as long as e > is sufficiently small. 

This concludes our discussion of the error matrix E(A). Wc only note two things at this point. Firstly, 
the bound H185() proves that 

(191) ipNix,t)=2i lim APi2(A) +0(7V-^), 

A — >oo 

as ^ oo because P(A) = E(A)P(A) for |A| sufficiently large, and h ^ Hn = A/N. Therefore, the strong 
asymptotics of the A-soliton arc provided by the modulated multiphase wavetrain that arises from the 
solution of the outer model problem for P(A) in terms of Riemann theta functions of genus G. In particular, 
the curves in the {x, i)-plane along which the genus changes abruptly are the caustic curves seen in Figures^" 
0] Secondly, we want to point out that the error estimate of 0{h) ~ 0{N~^) in H191|l is an improvement over 
the error bound obtained for the same problem in (10) . The improvement comes from (i) the e- modifications 
of the contours near A = which obviates the need for a local parametrix near the origin (this was also 
used to handle "transition points" in [J) and (ii) the careful tracing of the influence of the functions S{X) 
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and T(A) through the asymptotics, especially their explicit removal near the band/gap endpoints via the 
near-identity transformations P ^ Q. 



3.8. The formal continuum-limit problem. Much of the above analysis is based on the facts that 
5(A) = 1 + 0{h) and T{\) = 1 + 0{fi) under the assumptions in force on the relation between the contours 
on which these functions appear in the jump matrix and the contours C and Coo. These two functions 
measure the difference between the discreteness of the eigenvalue distribution and a natural continuum limit 
thereof (a weak limit of a sequence of sums of point masses). It is tempting to notice the role played by 
the approximations S{X) ~ 1 and r(A) w 1 in the rigorous analysis and propose in place of the problem 
for P(A) an ad hoc "continuum-limit" Riemann-Hilbert problem for a matrix P(A); the conditions of this 
Riemann-Hilbert problem are precisely the same as those of the Riemann-Hilbert problem governing P(A) 
except that in all cases one makes the substitutions S{X) 1 and T(A) 1. 

The rigorous analysis described earlier proves that, as long as the contour geometry admits the approxi- 
mations S{X) ~ 1 and T{X) w 1 one may also compute the asymptotic behavior of the A^-soliton by studying 
the formal continuum- limit Riemann-Hilbert problem for P(A). However, it turns out that there are some 
circumstances in which the conditions that constrain the contours of the Riemann-Hilbert problem for P(A) 
are inconsistent with the approximations S{X) w 1 and T(A) ~ 1. We will show that this is not merely 
a technical inconvenience standing in the way of analyzing the large N limit with the help of the formal 
continuum- limit Riemann-Hilbert problem for P(A), but that the modifications necessary to complete the 
analysis of P(A) rigorously under these circumstances introduce new mathematical features that ultimately 
provide the correct description of the secondary caustic. 

Note also that making the substitutions 5(A) 1 and r(A) 1 is completely analogous to the ad hoc 
substitution P{X) PiX), some of the consequences of which were described in S 13.2.21 These arguments 
suggest that extreme care must be taken in relating asymptotic properties of the matrix P(A) with those 
of the matrix P(A). Probably it is better to avoid analyzing P(A) and instead keep track of the errors by 
working (as we do in this paper) with P(A) directly. 

3.9. Dual interpolant modification necessary for contours passing through the branch cut. The 

nonlinear equations that determine the endpoints only involve quantities related to the continuum limit of 
the distribution of eigenvalues within the imaginary interval [0,Zj4]. Indeed, the rational function P(A) has 
been replaced by e^'-^^^^ at the cost of a factor 5(A), which is unifomly approximated by 1 on K-i. It 
turns out (see § that these equations admit relevant solutions that evolve in x and t in such a way that 
the contours become inconsistent with the fundamental assumption that the region D_i contains all of the 
discrete eigenvalues (poles of m(A) in the upper half-plane). This can happen because the function L(A) is 
defined relative to the contour C which while having the same endpoints as [0,«j4], is otherwise arbitrary, 
and the endpoint equations are analytic as long as the endpoint variables are distinct and avoid C. Thus, 
it can (and does) happen that a band / evolves in x and t so as to come into contact with the locus of 
accumulation of eigenvalues. When this occurs, we have to reconcile the facts that on the one hand the 
solution of the endpoint equations may be continued (by choice of C) in x and t so that the band / passes 
through the interval [0, iA] completely, while on the other hand the function 5(A) can no longer be controlled 
and the continuum-limit approximation can no longer be justified in the same way. 

The situation can be rectified in the following way. Returning to the matrix m(A) solving the discrete 
Riemann-Hilbert problem characterizing the A^-soliton, we remove the poles by taking into account three 
different interpolants of residues rather than just two. Consider the disjoint regions Di, D^i, and D^s in the 
upper half-plane shown in Figure |H1 We then remove the poles by going from m(A) to M(A) by exactly the 
same formulae as before, except in the region i?_3 and its complex conjugate image in the lower half-plane. 
In L)_3 we define instead. 



(192) M(A) m(A) 



for A e £)_3 , 



1 

-iP(A)e[2«Q(^)-3'^''(^)l/'^ 1 

and define M(A) for A G D*_^ by M(A) = cr2M(A*)*cr2. 

Once again it can be checked directly that M(A) is analytic at all of the poles of m(A) in the open domains 
_D_i and I?-3. The same argument works if it happens that a pole of m(A) lies exactly on the arc K^2 
separating the two domains. Next, we introduce the function .g(A), which we assume is analytic except on the 
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— e e 



Figure 8. The three regions Di, D-i, and and the boundary contour arcs Ki, Kq, 



contours K-i, Kq, and their complex conjugates, and we define a new unknown as before by setting 

N(A) := M(A)e~^*^^^°'^/''. On the arcs ATi, Kq, and A'_i the jump conditions relating the boundary values 
taken by N(A) are exactly as before, while on the arc Ar_3 we have 



(193) N+(A)=:N_(A) 
and on the arc A'_2, 

(194) N+(A)=N_(A) 



1 

iP(A)e[2»Q(^)-3»e''(A)-2g(A)]/?i I 

p-[s+(A)-9-(A)]/ri 



A e K-3 , 



iT{X)e 



[2jQ(A)-2,:e"(A)+L(A)-3+(A)-3_(A)]/ri g[s+(A)-s_ (A)]/n 



A e a:_ 



We choose the contour C, relative to which the functions A(A) and r(A) are defined, such that the contour 
K-2 forming the common boundary between the domains -D_i and D_3 lies entirely in the left half-plane 
according to C U Coo U C* U C^. For convenience we will take C to lie in the closure of Z)_i U Z)_3 and to 
connect A = to A = passing once through the first intersection point of A'_i and A'_2 (in the direction 
of their orientation, beginning with A = 0). We also assume (without loss of generality) that this intersection 
point turns out to lie in a gap of A'_i U Ar_2. See also Figure]^ 

The next step is to remove the jump discontinuities on the real intervals (— e, 0) and (0, e) by exactly the 
same factorization as was used previously (see (|81|) and the subsequent definition of 0(A) in terms of N(A) 
and the discussion thereof). Here the contour arc AT^ lies in the region D^i as before. 

The jump relations satisfied by 0(A) for A on Kq and AT-i are of exactly the same form as before (see 
(|91|l and (|92|l '). Furthermore, because A'_2 is in the left half-plane according to C U Coo U C* U C^, we may 
write the jump relation for 0(A) for A on A'_2 as 

e-[ff+(A)-s-(A)]/R 

iT(A)e[2*<3(^) + -^'(^)-»«"(A)-ff+(A)-3-(A)]/ri g[5+(A)--g-(A)]/n 



(195) 0+(A) = 0_(A) 



A e Kn 



We note here that the exponents of the jump matrix elements are of exactly the same form as for A G A'_i. 
Since T(A) = 14- 0(h) holds uniformly on A'_2 as long as we arrange that this contour is bounded away from 
the endpoints of C, we may asymptotically analyze the Riemann-Hilbert problem for 0(A) by choosing .g(A) 
in exactly the same way as we did earlier, with the same integral conditions on the endpoints of the bands. 
In particular, the characteristic velocities in the Whitham equations arc exactly the same analytic functions 
of the elements of v as in the simpler configuration. The contours of discontinuity for P(A) in the modified 
configuration are shown in black in Figure El 

33 




Figure 9. The discontinuity contours for P(A) in the modified configuration. The arc 
connecting Ai and A2 between dflf is the band Ii which may now cross the interval [0, iA] . 

With regard to the determmation of the function g{X), the only essential difference between this approach 
and the former approach is that here a new inequality is required to hold for A S in order that the 
corresponding jump matrix for 0(A) be a small perturbation of the identity matrix. From H193() . we sec that 



(196) 



0+(A) = 0_(A) 



A G if_3 . 



1 

i5'(A)e[2»Q(^)+-^(^)-»»"(^)"2s(A)-2ie''(A)]/;i I 

The exponent is an analytic function in I?-3, and taking a boundary value on /i_2 we see that 

2iQ{X) + L(A) - ie°{X) - 2g{\) - 2ie°{X) = 2iQ{X) + L{X) - ie°{X) - 2.g_(A) - 2i9°{X) 

(197) 

= (/)(A) - i6i(A) - 2i6i"(A) . 

By assumption, the arc K^2 contains a part of a gap, so we may evaluate the boundary value in that gap, 
where 0{X) is a real constant and find that the relevant inequality is 

(198) 5R(2iQ(A) + L(A) - t6i"(A) - 2.g(A) - 2i0°{X)) = 5R(0(A) - 2i6i"(A)) < , AG A'_3 . 
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Here (f>{X) denotes the analytic continuation from any gap on K-2 to (while different analytic functions 
for each gap on K-2, they all differ by imaginary constants which play no role in H198() ). 

The term 23(0° (A)) represents a nontrivial modification to the inequality 5R(0(A)) < that must be 
satisfied in the gaps of A"_i U K^2- As wc will see in § |S1 the violation of this inequality is exactly the 
mechanism for the phase transition from G — 2 (genus two) to G = 4 (genus four). It should be stressed 
that this inequality is fundamentally an artifact of the discrete spectral nature of the A^-soliton. In other 
words, if we simply propose (as described in 13. 8|) a model Riemann-Hilbert problem for an unknown P(A) 
in the original scheme with all jump matrices modified simply by setting S'(A) = T(A) = 1, we (i) will not 
care about the implications of a portion of the contour K-i meeting the interval [0,«A] (which has meaning 
only as an artifical branch cut of an analytic function) and (ii) will only need the inequality 3fi(0(A)) < 
and ascribe no particular meaning to the modified quantity 3?((/)(A) - 2i0°{X)). 

4. The Genus-Zero Region 

4.1. Validity of the genus-zero ansatz. Ansatz failure and primary caustic. In Chapter 6 of 
the genus-zero ansatz for the g-function of the iV-soliton is constructed and it is shown that the inequality 
3?((/)(A)) < is satisfied throughout the sole (terminal) gap for < t < ti{x), where the curve t = ti{x) 
is the primary caustic curve asymptotically separating the smooth and oscillatory regions in Figures 
The primary caustic curve is described by obtaining the endpoints Ao(x,i) and XQ{x,t) from the genus-zero 
ansatz for .g(A) and then eliminating A from the equations 

(199) 4>'{X)^0, 3?(</)(A)) = 0. 

In other words, the first equation says that A ~ X{x,t) is a critical point of 4>{X) with endpoints Xa{x,t) 
and XQ{x,t), and then the second equation is a real relation between x and t, that is, a curve in the {x,t)- 
plane. The equations 1199|) describe the existence of a critical point A for (j){X) that lies on the level curve 
3?(0(A)) = 0, which is the boundary of the region of existence for the terminal gap Fi. The existence of such a 
point on the level indicates a singularity of the level curve, generically a simple crossing of two perpendicular 
branches. In other words, if (5 > is small then for ti{x) — 5 < t < ti{x) the region where 5R(^(A)) < holds 
is connected albeit via a thin channel, delineated by approximate branches of a hyperbola, through which 
the gap contour Fi must pass, whereas for ti{x) < t < ti{x) + S the region where 3?((/)(A)) < holds becomes 
disconnected, and it is no longer possible to choose Fi so that 3?(^(A)) < holds throughout. When t = ti{x) 
the hyperbolic branches degenerate to lines crossing at A = X{x,t). This mechanism for the formation of 
the primary caustic is illustrated in the spectral A-plane in Figures 6.11 and 6.12 of |1()) . In particular, it 
is known in the case of the A^-soliton that ti{x) is an even function and inix£Rti{x) = ii(0) = {2A)~^. 
In ^Oi it is shown that if \x\ decreases with t fixed from a point (xcrit, ^crit) at which H199() holds, then the 
genus-two ansatz takes over. A new band is born with two new endpoints Xi{x, t) and A2(a;, t) emerging from 
the double point A where the singularity of the level curve for 5ft((/)(A)) = occurs. Therefore, the primary 
caustic is a phase transition between genus zero and genus two. 

In fact, the conditions H199|) are not specific to the failure of the genus-zero ansatz, but may be considered 
for any genus, being as they simply describe the change in connectivity of the region where the essential 
inequality 5R(0(A)) < holds. Thus, we may consider the conditions p99|l for general even genus with 
the aim of determining whether phase transitions beyond the primary caustic can occur for the iV-soliton. 
Simultaneous solutions of the two equations in 1)199(1 can be interpreted as collisions between branches of the 
level curve 5R(0(A)) = 0, and these branches may emanate from band/gap endpoints Afe(x,t), or from the 
branch cut C, or they may be noncompact. (As level curves of a harmonic function they cannot close on 
themselves unless they enclose a singularity of 4>{X).) As shown in Figures 6.11 and 6.12 of ^Hli the primary 
caustic is caused by the collision of a noncompact branch coming in from infinity with a branch joining 
Xo{x,t) to the origin. It is not difficult to analyze the function 0(A) for general even genus to determine its 
singular points and asymptotic behavior as A — > oo, and such analysis leads to the conclusion that for genus 
two no further branches exist for large A that can eventually collide with branches connecting band/gap 
endpoints. This fact suggests that further phase transitions (higher-order caustics), if they exist, may come 
about for a different reason than the primary caustic. The main point of this paper is to show that a 
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secondary caustic indeed occurs for the iV-soliton in the large N hmit, and that its mechanism is indeed 
different from that leading to the formation of the primary caustic. 

The validity of the genus-zero ansatz for g{X) for small time is in fact established in ^U] for a general 
class of "reflcctionless" initial data for the focusing nonlinear Schrodinger equation that leads to a discrete 
Riemann-Hilbert problem for m(A) of the sort relevant in the study of the iV-soliton. Similar results have 
been found by Tovbis, Venakides, and Zhou | 15| for certain special-function initial conditions that are not 
reflcctionless. 

4.2. Fourier power spectrum and supercontinuum frequency generation. In applications of the 
focusing nonlinear Schrodinger equation to nonlinear fiber optics, the limit 4 is interesting because it 
corresponds to weak dispersion, a situation in which nonlinear effects arc dominant over linear dispersive 
effects. While conservative linear effects leave the Fourier power spectrum of a signal unchanged, strongly 
nonlinear effects can change the spectrum signiflcantly, and one application of weak dispersion in fiber 
optics is the generation of broadband spectra from input signals with a narrow power spectrum. This is 
known in the literature as supercontinuum generation. While it is sometimes thought that supercontinuum 
generation occurs as a result of the wave breaking at the primary caustic, in fact the effect is already present 
for t < 'mfr^^]^ti{x). Indeed; in this region, the A^-soliton has, for each x G R the representation Ijiyifl 
corresponding to the genus-zero ansatz for .g(A). Moreover (see jl(J| ') in the genus-zero case, 

(200) 2i lim APi2(A) h^e-"'"''' , 

\ — >oo 

where Ao = ao{x^t) + i6o(x,i) and hiQ — n^ix^t) are related by 

(201) ^ = 2„„. 

For the A^-soliton with < > 0, the endpoint coordinates ao(a;, t) and 60(2;, t) satisfy 

J2 _ u2 ^ ^-2ui 

(202) al 



/12 _ / 2 , J-2L4 



and 

(203) X = -2taQ + (^arcsinh (^ ""^'^ ^ ^ . 

Note that if < t < {2A)~^ = inix^R ti{x), the relation 12U2|) defines a curve in the real (ao, 5o)-plane with a 
compact component having a "figure-8" shape with two lobes joined at the origin, one enclosing the interval 
(0, iA) and the other enclosing the interval {—iA, 0). For such t, the correct solution (in the sense of leading 
to topological conditions for contours in the complex plane admitting the genus-zero ansatz) of (|202|) and 
(|203|l necessarily lies on the upper lobe of this compact component, on which it is uniquely determined by 
(|203|l as a function of x. In other words, (|203|) parametrizes the upper lobe by x € R, taking M_ to the left 
half of the lobe, to the right half, and a; = to the point above iA where the upper lobe intersects the 
imaginary axis. As x ^ —00, Ao(a;,t) tends to zero from the second quadrant, and as a; ^ -|-c», Xo{x,t) 
tends to zero from the first quadrant. See Figure [TUi 

The power spectrum of tpNiXjt) in the globally genus-zero regime (for t < (2^)'^) can be approximated 
by using the method of stationary phase to analyze the Fourier transform of the approximate solution (|200|1 . 
To use the method of stationary phase to study the approximate power spectrum 

-1 />oo 2 

— / 6o(x,t)e-'''=°(^'*)/''e-*"^dx 



(204) F{uj,t) 
it is useful to introduce the scaling uj = il/h and consider fl fixed as ft J, 0. Then, 

(205) F{uj,t) 



1 r 

and the stationary phase points satisfy 

(206) =^ + n^2ao(x,t) + n. 

ox 
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Figure 10. The upper lobe of the compact part of the curve (f^n^ shown for At = 
0, 1/8, 1/4, 3/8, 1/2 /rom inside out. Fort > the relation H2()3(l parametrizes the curve 
from the origin to itself in the clockwise direction. 



According to the parametrization of the upper lobe of the compact component of the curve H202() given by 
()203|l . there exists a finite number M > (the maximum of 2|ao| over the upper lobe) such that 

• There are no stationary phase points for |ri| > M. 

• There are two simple stationary phase points x+{fl) < X-{^1) < for < < M. 

• There are two simple stationary phase points < X-{^1) < x+{^l) for —A/ < < 0. 

• For = M there is a double (degenerate) stationary phase point. 

• For fi = there is only one finite simple stationary phase point at a: = but the phase is also 
stationary at a; = oo. 

Thus, the power spectrum is negligible for frequencies uj of magnitude greater than M/h. The scaled cutoff 
frequency M depends on both A and t, but from dimensional analysis of H2U2|) it can be seen that M/A is a 
function of the combination At alone. The scaled cutoff frequency is plotted in Figure ITTI over the interval 
< At < 1/2 in which the genus-zero ansatz holds for all a; G R. 




Figure 11. The scaled cutoff frequency is a monotone increasing function of At. 



In the bulk of the spectrum |f2| < M, the power spectrum is proportional to h. Indeed, the relevant 
stationary phase formula for the generic case of |f2| < M and for nonzero fixed 17 is 



(207) 

where 

(208) 



F{u;,t) 



h 

An 



A+{nf + A^{nf - 2A+{n)A^{n)sm f ) 



V ^ / 



0(ft2), 



A±in) :-- 



bo{x±in),t) 

y/\d^aoix±{n),t)\ 
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and 



(209) := -KQ{x+{n),t) + Ko{x^{n),t) - nx+{n) + nx-{n) . 

Considered as a function of Q in the bulk \^\ < M, the power spectrum is a rapidly oscillatory function such 
that 

(210) ^iA+{n)-A_in)f + oite)<F{uj,t) < ^{A+{n) + A.{n))^ + o{tf). 

The rapid oscillation in the spectrum is a kind of "ripple" that has a characteristic spectral width on the 
scale given by 

li h 

(211) An^ — . 

^ ' \x+{n)-x^{a)\ 

Since |a;+(ri) — x_(r2)| tends to zero as |r2| M and grows without bound as ^ 0, the ripple becomes 
less pronounced near the edge of the spectrum and more violent near the center of the spectrum. Near the 
edge of the spectrum the characteristic spectral width still tends to zero with h but at a slower rate; in this 
regime more detailed stationary phase analysis of the power spectrum would reveal Airy-like behavior on a 
scale where |ri| — Af is proportional to an appropriate fractional power of h. In any case, uniformly over the 
region < M, the ripple is rapid enough to give meaning to the local average (weak limit with Q fixed): 

(212) {h-^F{uj,t)) := [A+{n)^ + A-{n)^] . 
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The weak convergence of the power spectrum to a broad plateau of spectral width (on the uj scale) of h^^M 
can be considered as evidence of supercontinuuni frequency generation for the A'^-soliton in the genus-zero re- 
gion, before any wave breaking occurs. Furthermore, these calculations show that the broad spectral plateau 
develops in absence of the Raman effect or higher-order dispersion, which are not part of the basic nonlinear 
Schrodinger model but are frequently thought to play an important role in supercontinuuni generation. 

5. Numerical Computation of Endpoints and Contours in the Genus-Two Region 

Our main interest in this paper is to analyze the A'^-soliton in the genus-two region to determine and 
explain the mechanism for the secondary caustic curve that can be observed in Figures In the genus- 
two region, there are three complex endpoints Ao(a;, t), Ai(a;, t), and A2(a;, t) that satisfy the six real equations 

Vi = 

i?i = 

Mo = 

Ml = 2t 



(213) 



Ma = x + 2t^ 



2 

a„ 

n=0 



2 2 ^ ^ s 2 2 

M3 = a;^a„ + 2t^ --62j +2i^ ^ a™a„ . 

n— n—0 n— Om— n-t-1 

We wrote a Matlab code that continues a solution of these equations along a path from point (a;o,to) at 
which a solution is known for which the Jacobian of the system is nonsingular to a nearby point {xi,ti). 
To have a starting point, we need to choose a point (xcrit, ^crit) on the primary caustic, to which there 
corresponds a single endpoint Ao(xcrit, ^crit) of the genus-zero ansatz, and a critical point A at which H199() 
holds. In |10l it is shown that for x = Xcrit and t = t^rit the equations (|213|) admit the solution given by 
A2 = Ai = A and Aq taking the same value as in the genus-zero case. Since the Jacobian of H213() is singular 
when endpoints coalesce, we need to begin with an analytical perturbation calculation to move slightly away 

38 



from the primary caustic into the gemis-two region. We did the calculations for Xci-a = and tcrit = (2^) ^, 
where Ao = A = and analytically determined asymptotic formulae for Ao(a;,t), \i{x,t), and \2{x,t) in 
the limit t | {2A)~-^ with a; = 0. With the help of these formulae, we could employ our continuation code 
beginning with a point slightly beyond the primary caustic where the three endpoints are distinct. Thus, 
we may obtain numerically, for any (x, t) beyond the primary caustic, the analytic solution of the equations 
()213|1 that correctly emerges from the primary caustic. 

We also wrote a Matlab code that takes a point (x, t) beyond the primary caustic and the corresponding 
endpoints Ao(a;, t), Ai(x, t), and A2(x, t) as input, and computes the band /q connecting Ao(x, t) to the origin 
and the band Ii connecting Ai(a;, t) and A2(x, t) by integrating numerically the equation ^{R{rf}Y{rj) drf} — 
to determine the level curves. The code also finds the regions of the complex plane where the inequality 
5R(0(A)) < holds; the gaps Fi connecting Ao(x, t) to Ai(x, t) and r2 connecting A2(x, t) to an e neighborhood 
of the origin must lie in these regions. By plotting this information, we can visualize the dynamics of the 
genus-two ansatz in the spectral plane, and hopefully determine the nature of the secondary caustic. 

A representative sequence of such plots is given in Figure E| Here we see the bands of the G = 2 contour 
in red and the region where the inequality 3?((/)(A)) < holds shaded in blue. For the contour C we simply 
take the imaginary interval [0,i^], shown with a black line. Each spectral plot is paired with a plot like 
Figure^lbut with a small red dot indicating the corresponding point in the {x, <)-plane. The first plots show 
the genus-two configuration for a point in the (x, t)-plane just beyond the primary caustic. The band Ii in 
the left half-plane has just been born from a point A where ()199|l holds. The bands can clearly be joined 
together with gaps Fi and F2 lying in the shaded regions to complete a contour Kq U K-i that encircles 
the locus of poles for m(A). Notice, however, that according to the subsequent plots in this figure the band 
Ii approaches the locus of accumulation of poles for m(A) as t increases. Although a collision of one of the 
endpoints with the interval [0, iA\ would imply a failure of the genus- two ansatz in its original formulation 
(that is, without the modification described in § I3.9|l . it is clear by comparison with the corresponding 
(x, t)-plane plots that the secondary caustic is still far off. 

Figure IT!^ continues the evolution. Note that between the first two frames, the band Ii collides with the 
interval [0,i^], so we must switch to the modified version of the genus-two ansatz described in 5 13.91 Thus, 
a new contour arc A'_3 must be included and a new inequality 5R(0(A) — 2*6'" (A)) < must be satisfied 
along this contour. The region where this inequality is satisfied is shaded in green. Thus, the contour arc 
K-^ must lie in the green shaded region, while the gaps of K-i U K-2 must lie in the blue shaded regions. 
Furthermore, the contour C must from this time onward be deformed somewhat to the right to admit the 
passage of the band endpoint Ai(x,<) "through the branch cut". The deformed contour C we choose here, 
somewhat arbitrarily, is the union of the two indicated solid black line segments. It is clear that in each case 
the bands may be complemented with gaps in the blue regions to form (along with a small interval (— e, 0)) 
a closed loop. In cases where this loop necessarily cannot enclose the poles of m(A) accumulating in [0, iA\, 
it is clear that a curve iiT-s lying entirely in the green shaded region may be added that passes over the top 
of [0, iA\ and forms along with parts of the above closed loop the boundary of a region containing the poles. 

In Figure ini we continue the evolution toward a curve in the (x, t)-plane at which the asymptotic behavior 
abruptly changes for the second time, where we expect to find a mathematical explanation for the secondary 
caustic. Note that between the last two snapshots there is a pinch-off of the green shaded region where 
3?(0(A) — 2i0°(A)) < 0. When this happens, a contour lobe K-^ that serves to bound the region containing 
the poles of m(A) (soliton eigenvalues) and also along which the inequality 3?(0(A) — 2i6'°(A)) < holds can 
no longer exist. This is the mathematical mechanism for the development of the secondary caustic. 

To demonstrate more convincingly that the secondary caustic in the numerical reconstructions of the 
A'^-soliton solution for large TV actually corresponds to the curve in the (x, t)-plane along which the region 
in which the inequality ^{(j>{X) — 2i6^{\)) < is pinched off at a point, we used a secant method to find, for 
a given x, the i- value at which 5R((?!)(A) - 2ie°{X)) = where A is a complex root of (^'(A) - 216*°' (A). Level 
curves of 5R(0(A) — 2i6'°(A)) passing through simple critical points like A necessarily exhibit a characteristic 
perpendicular crossing singularity. This method yields (approximate) coordinates of points in the (x, t\ 
plane at which the configuration of level curves makes a transition between the two types of configurations 
illustrated in the final two frames of Figure ^1 These points lie along a curve t ~ t2 (x) that defines the 
location of the secondary caustic. Configurations corresponding to three of the points we obtained in this 
way are illustrated in Figure IT^ 
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Figure 12. Evolution of the contour configiffation for G = 2. Here x = 0.5 and t = 0.7 
through t = 0.9 in steps of At = 0.1. 5*66 the text for a full explanation. 
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Figure 15. Configurations in which the regmn where the inequality 3?(0(A) — 2i9^{X)) < 
holds is "pinched- off " delineate the secondary caustic. 



6. Transition to Higher Genus 



6.1. Violation of the new inequality. According to the numerics explained in § |31 the failure of the 
genus- two ansatz at the secondary caustic corresponds to the equations (compare with H199|) ): 

(214) (/)'(A)-2i6i°'(A) = 0, di{(f>{X)-2ie°{X))^Q. 

We expect the critical point A arising from these equations to be the origin of a new band, and to show this 
we need to develop the genus- four ansatz within the framework of the modified approach given in § 13.91 but 
where the new band appears on the contour This gives the equations implicitly defining the endpoints 

as functions of x and t a different appearance than if the new band were assumed to lie on K-i U if_2. 
These equations turn out to define a previously unknown branch of solutions of the quasilinear system of 
Whitham equations (as defined in S IH.4.1|I . 



6.2. On and beyond the secondary caustic. 

6.2.1. Endpoint equations for genus four with a band on The modification described in S 13.91 leads to 

a new inequality that must be satisfied for A € K-^. Moreover, the numerical computations in §|S1 indicate 
that the failure of the genus- two ansatz is due to the violation of this new inequality. We then expect, 
similarly as in the case of the transition from genus zero to genus two across the primary caustic, that the 
pinching off of the inequality region 5R(^(A) — 2i6'''(A)) < leads to the birth of a new band and a phase 
transition. 

In this section we describe the modifications to the endpoint equations that arise when we allow a band 
to lie on the arc K-^. In particular, we suppose that there is a single band in each of the arcs Kq, K-2, and 
ii'_3, and we develop equations corresponding to (|131|l and H133|l - (|134|) in this special case. The resulting 
contours of discontinuity for P(A) in this situation are illustrated in Figure [TCI To begin, consider the 
analytic matrix M(A) obtained as in S 13.91 from the triple interpolation procedure. As before we introduce 
the function g{X), but now we suppose that g{X) is analytic except on A'o, if-i, and K-^ along with 

the complex conjugates of these arcs. We then define a new unknown via N(A) := M(A)e-9(^)'^3A_ All 
the jump conditions relating boundary values are as before with the exception of the condition regarding 
boundary values on A'_3. In place of H193|) . wc find 

r e.-^9+W-a-m/h Q 

(215) N+(A)=N_(A) .p^y^^[2tQ(X)-Zie\\)-g^(\)-g^{\)]/h g[3+(A)-g_(A)]/ri ■ 

Continuing to follow the steps of our earlier analysis, we arrive at a Riemann-Hilbert Problem for an unknown 
matrix 0(A) which has no jump discontinuities on the real intervals (— e, 0) and (0, e), satisfies the jump 
conditions l|91() - H92|l . H195|l . and also satisfies 

iS'(A)e[2»Q(^)+-f^(^)-»«°(^)-f+(-^)-9-(^)^2je°(A)]/?> g[s+(A)-3_(A)]/rx 

From (|216(l . the band and gap conditions for the arc Ar_3 arc slightly different from those in H96|l and (|97|l 
which continue to hold in the arcs A'o, A'_i, and K-2- Namely, for a band in Ar_3, we require 

(217) 0(A) — 2iO^{X) = imaginary constant, 0{X) ~ real decreasing function, 
while for a gap in A'_3, there must hold 

(218) 3?(0(A) - 2i6l°(A)) < 0, e{X) = real constant. 

Then, we note that ^'(A) must satisfy H98() in the bands not on A'_3, (|99|l in all gaps, (|l()()(l . (|1()HI . and 

(219) g'+{X) + g'_{X) = 2iQ'{X) + L'{X) - 3ie°'{X), for A in the band on A'^3- 

These conditions amount to a scalar Riemann-Hilbert Problem for g'{X). Wc (re)define C/ to be the contour 
Ko U A'_i U K_2 U (-e, 0) U C, and we define /(A) as before by ifT?!^ and JTHSll- We then make the change 
of variables 

(220) h{X):=g'{X) + f{X). 

We see that the function h{X) is characterized by 1)1 13|) for bands not in -A-3, (I114)l - I|116() . and 

(221) h+{X) + h-{X) ^2iQ'{X)~3ie^"{X), in the band on 3. 



(216) 0+(A)=:=0_(A) 



for A e A'_3. 




Figure 16. The discontinuity contours for P (X) in the genus-four configuration with a band 
on the outer contour lobe K^^. 

Wc have thus traded our scalar Riemann-Hilbert Problem to determine the unknown g'{\) for one in which 
the unknown is h{\). We suppose that the endpoints of the bands 

Aq, Ai, . . . , A4 

are known, and we define R{X) by 

4 

(222) i?(A)2 = I](A-A„)(A-A:J, 

ri=0 

and the conditions (i) the branch cuts are the bands in A'o, K^2, and A'_3 and (ii) i?(A) = A'"' + O(A^) as 
A — > 00. We denote the bands by /q, Ii, and l2- Note that I2 has the same orientation as if_3. Finally, we 
are able to solve for h{\) by writing it as the product of R{\) and a new unknown function fc(A) as in (|119|l . 
From the definition of i?(A), the conditions on h{X) translate to conditions on fc(A) as in (|120|I - H124|I with 

(223) fc+(A) - fc-(A) ^ ^"3'^^) for A in the band on /v_3. 

R+{\) 
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Now, using the differences of boundary values of fc(A) in (|120|l - p22|l and (|223|l . we can write down an explicit 
formula for fc(A). We define fcc/(A) := klj\x) + fc[f''(A) where 



(225) k\f\X):^^- 



2i 7c/\bands Riv)iV~ 



We note that we can rewrite fc^,^-* as 



)'(r;) _^ 3 /■ 



l!7(A) 

Finally. fc(A) must have the following form: 

(226) k(X):=ku{X)-ku{X*)*. 
A residue calculation shows that 

fcW(A) - kij\xr = '^+iu{x)-iu{xr, 

hence 

(227) fc(A) = +IuiX)~IuiXr + kl^\x) - k'^\xr. 

Finally, we derive the analogues of by enforcing the decay condition (|124|1 . We expand the integrals 

appearing in (|227|l for large A. We note that 6*"' (A) = in, so 



-in) 



drj 



\-^-P + 0{X-'^) as A ^00. 



Starting from the notation (|129|1 . wc introduce modified moments 



(228) Mp(Ao, . . . , A4) := A/p(Ao, . . . , A4) - 35R ^ drjj , 
and we thus find that 

(229) = Mp(Ao,...,A4), for p = 0,1, 2, 

(230) 2t = M3(Ao,...,A4), 

4 

(231) a; + 2i^a„ = M4(Ao,...,A4), 

n=0 

(232) a;^a„ + 2t ('a2_i52V ^ ^ a„a„ J = M5(Ao, A4) . 

n=0 \n=0 ^ ^ n=0 m=n+l / 

These equations amount to real constraints that must be satisfied by the 5 complex endpoints in the upper 
half plane. 

To ensure that 5(0(A)) = in the gaps, we require as before that (recall that An is a small counter- 
clockwise oriented closed contour surrounding the band /„) 

(233) 5R^^ g'(7?)dry^ = 0, n=l,2. 
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Additionally, we need to enforce that 3?((^(A)) = for A in the bands Iq and /i and that di{4'{X) — 2i9'^ {X}) ~ 
for X G l2- The condition 



□0// 



(234) 3? (^y^ [2g'(77) - 2zQ'(ry) - i'(r?) + lO^'ir^)] drj j ^ , 

where the contour of integration is arbitrary in the domain of analyticity of the integrand, will guarantee 
that (/)(A) is purely imaginary in the band Ii. On the other hand, we know that (/)(A) — 2i9^(X) is guaranteed 
to be constant on I2 by our construction. To ensure that this constant has zero real part, we require that the 
net change of the real part of ^(A) across the second gap is equal to the real part of 216*" (A3). The second 
vanishing condition is 

(235) 3? ^ ' [2g'{r,) - 2iQ'ir^) ~ L'ir,) + tO^'ir,)] dr, - 2*0"(A3)^ = , 

where again the path of integration is arbitrary within the domain of analyticity of the integrand. Together 
with H229() - (|232|l . equations (|233|I - H235() make up 10 real equations for the real and imaginary parts of the 
endpoints Aq, . . . , A4. 

The conditions can be written in a simpler form that is useful for subsequent calculations. 

Define 

and then set Y{X) := ^c/(A) — Yij{X*)* . Therefore Y{X) is a function that is analytic where Y{X) is with the 
exception of the contours A2 and Aj. Then, the conditions H233() arc equivalent to 

(237) i?„(Ao, . . . , A4) 3? ^ Riv)Yii^) d^i^ = , 
for n = 1 and n ~ 2. Similarly, the condition 1234|) is equivalent to 

(238) l7i(Ao,...,A4) :-3? ( / ' i?(ry)f (t?) d?/) =0. 



'A 



2n~2 



In all of these formulae, the contour of integration is arbitrary as long as it lies within the domain of 
analyticity of the integrand. Let Am denote a point on A2 . Then, the condition l)235|l may be written as 



/ R{Tj)Y{r])dTj+ R{Tj)Y{rj)dT^ + 2i9°{X2)-2Tr ^7^ =0. 

JX2 JXm Jx2 J 

Here the paths of integration lie in each case within the domain of analyticity of the integrand. 



6.2.2. Birth of a new band upon crossing the secondary caustic. Here we show that when (a;, t) is a point 
on the secondary caustic, so that the equations (|214|l hold, the equations derived above that determine the 
endpoints of the genus-four ansatz with band I2 on the outer contour lobe A'-a admit the solution where 
Ao(a;,t), Xi{x,t), and X2{x,t) are all equal to the corresponding endpoints for the genus-two ansatz at the 
same point in the (a;, i)-plane, while X3{x,t) = Xi{x,t) = X. 

To show this we first note that the conditions defining the function g'{X) in terms of the endpoints yield 
the same function in both the genus- two case (ignoring A3 and A4) and in the genus- four case described above 
(assuming A3 = A4). Since Aq, Ai, and A2 arc such that g'{X) = 0{X~^) as A ^ 00 (because the genus-two 
version of 5' (A) satisfies the appropriate suite of four moment conditions), we also know that the genus-four 
version of g'{X) satisfies the corresponding suite of six moment conditions (see (|229|l - (|232(l 'l. 

It therefore remains to verify the remaining integral conditions H233|l - (|235|) involving g'{ri). Since the 
function g'{ri) is the same whether we consider the genus- two or degenerate genus- four version, the conditions 
Ri — Q and Vi = arc automatically satisfied, as these are equivalent to the conditions Ri = and Vi = 
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satisfied by the genus-two version of .g'(?7). Furtliermore, since g'(ri) is analytic in a neighborhood of A, we 
have by Cauchy's Theorem 



(240) f g'{Tj)dv = 0, 

which in particular implies that the real part is zero, so i?2 = is satisfied as well. Finally, we observe that 
the second vanishing condition H235|l is, when the point A is taken as the second endpoint of the gap 
satisfied as a consequence of the definition H214|l of the critical point characteristic of the secondary caustic. 
This completes the proof that each point on the secondary caustic, as defined by (|214|) for the genus-two 
ansatz, corresponds to a degeneration of an admissible genus-four configuration with A3 = A4 = A. 

Showing that the points A3 and A4 separate in a direction admitting a new band I2 between them as t 
increases beyond the secondary caustic is slightly complicated by the fact that the Jacobian of the system of 
endpoint equations is singular in this case. However, a direct calculation shows that the 8x8 Jacobian of all 
endpoint equations with the exception of (|233|l for n = 2 and (|235|l taken with respect to the independent 
variables Aq, Ai, A2, (A3-|-A4)/2, and their complex conjugates is nonsingular in the degenerate configuration 
when the remaining independent variables S := (A4 — A3)/2 and S* :~ (A| — A3)/2 both vanish. This 
reduces the problem to the study of two equations and two unknowns. It follows that near the degenerate 
configuration, the eliminated variables are all analytic functions of x, t, S, and S*, and it can be checked that 
all are invariant under 5 —S and 5* —S* independently. This actually makes them analytic functions 
of X, t, 5^, and 5*'^. 

The remaining equations are i?2 = and V2 = 0, and the left-hand sides of these are now to be considered 
functions of x, t, S, and S* . Knowing that they admit the solution 6 = 5* ~ when {x,t) is a point on the 
secondary caustic, we wish to solve them for the variables S and S* for nearby {x,t). We will show that there 
is a unique solution satisfying the condition that S and S* are complex conjugates of each other for real x 
and t (there are other solutions that bifurcate from 5 ~ 5* ^ that do not have this symmetry, which is the 
primary source of difficulty). We define polar coordinates by 

(241) ,5 = re'^ , 6* = re"'^ , 

and consider the equations R2 = and V2 = naturally complexified (that is, the left-hand sides extended 
as analytic functions of four complex variables x, t, r, and 9. Being as the complexification of R2 can be 
written as the average of two integrals, one on a path surrounding the two endpoints A3 and A4 and the 
other on a path surrounding the two endpoints A3 and A|, it is easy to see that R2 is an even function of S 
and S* individually that vanishes when 6 = 6* =0, and therefore the equation R2 = has the form 

(242) = a{x, t)r'^ cos(2(0 — (3[x, t))) + quartic (and higher order) terms in r. 

This equation clearly defines two analytic curves crossing orthogonally at the origin in the real polar plane, 
with angles 9 = (3{x, t) + 7r/4 + mi/2. 

Now, the complexified form of V2 defines a multivalued function of 6 and 6* , even when we restrict to the 
real subspace where these two variables are complex conjugates of each other. This comes about because the 
value of the integral terminating at A3 in V2 can change its value if A3 and A4 are permuted. This monodromy 
effect can, however, be removed in the real subspace, because the deformation (<5, 5*) {6e^'^ ,6*e~^'^) for 
T going from zero to tt (which exchanges the endpoints and their conjugates) simply adds to V2 a copy of 
the function i?2- This suggests that in place of V2 = we can solve an equivalent modified equation that is 
single-valued in the real subspace. Indeed, we opt to replace V2 = by vf" = where 

Vf := Riv)Yirj) dr^^ - ^ log (^A^ • Riv)Y+iv) dr^ 

(243) 

+ / Riv)Y{ri) drj + 2i9^{\2) - 2tt dry . 



Here Y+{t]) refers to the boimdary value taken by Y{ri) on A2 from the interior of this closed contour. 
Clearly the equation = has the same solutions as does V2 = 0, subject to the condition R2 = 0. But 
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the function Vf has the advantage of being single- valued in the real polar plane of (r, 9). Moreover, we may 
;e vf in the form 

V* = ^ij^ ' R{v)Y{v) drj - ^ log(J) £ R{v)Y+ (,7) drA + ^ log(r) • 3 Ril^)Y+ (r;) dri^ 



(244) 



/ R{ri)Y{T])dr] + 2i9°{X2)-2Tr dr] 
V Aq A2 > 



Here we see that the first term of Vf is complexified as an analytic function of S'^ and 5*"^ individually. 
Therefore it has an expression in the real polar plane of the form f{x, t) + O(r^). The complexification of the 
second term, while single- valued when 6 and S* are complex conjugates of one another, is not single- valued in 
a full neighborhood of (5 = (5* = 0. Therefore it cannot have a two- variable Taylor expansion; however the 
factor multiplying log(r) clearly is an even analytic function of both S and S* that vanishes when S = S* = 
(this is of course just the complexification of the harmonic conjugate i?2). Consequently, in the real polar 
plane we have 



(245) V* = ^^r^ log(r) sin(2(0 - P{x, t))) + 7(x, t) + O(r^) , 



where ^{x, t) is an analytic function of x and t that vanishes to first order on the secondary caustic; it is 
the value of given by H244|) evaluated with all endpoints in the degenerate configuration. Equating this 
to zero, we see that the balance must occur between the log(r) term and the constant term, so the the 
distance from the secondary caustic scales as log(r) as the new band is born. The curve described by 
(|245|l near the origin in the polar plane consists of two approximately hyperbolic branches with asymptotes 
rotated from the curves described by the equation i?2 = by ■n/A. Thus, we obtain for each point {x,t) 
sufhciently close to the secondary caustic a pair of opposite solutions for 5 = re'^. As ^{x,t) changes sign, 
the two points come together aX 5 — Q and then go back out along the perpendicular direction. 

Determination of the sign of 7(0;, t) analogous to what is carried out in |10| is then needed to show that 
the points X^^ix^t) and A4(a;,t) separate in the direction conducive to the completion of a genus- four ansatz 
as long as x decreases (or t increases) as the secondary caustic is crossed. Combining these steps completes 
the proof that upon crossing the secondary caustic the asymptotic description of the iV-soliton is given in 
terms of Riemann theta functions of genus four. 

We also note that the equations i?2 = and V2 = Q imply that the new band I2 and the boundary 
curves of the region where the inequality 3?((/)(A) — 2i6^{X)) < holds are trajectories of the quadratic 
differential i?(?7)^yi|_(77)^ dr/^, while the other two bands and the boundaries of the region where the inequality 
9fi((?!)(A)) < holds are trajectories of R{r])^Y^{7])^ drf. 



6.2.3. Whitham equations. Recall from i) l3.4.1l that for arbitrary even genus G (including G = 4 as a special 
case), if the endpoints occur along A'o U K^i or Kq U U K^2, then the endpoints implicitly defined 
as functions of {x, t) via the relations Ei = ■ ■ ■ ~ -E2G+2 = satisfy a universal set of quasilinear partial 
differential equations, namely the Whitham equations H159() with complex characteristic velocities given as 
functions of the endpoints only via (|160(l . 
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We have just seen, however, that with 
El := Vi, 



(246) 



E2 
E3 

En 

Eg 



i?l , 
V2, 
i?2 , 
Mn-5 
M3- 



for n = 5, 6, 7 , 



2t, 



Eq 



E- 



10 



Ma 



2t a„ , 



n=0 



4 / 4 . 1 \ 4 

:^a„ - 2t I ^ f - -62 j +^ ^ a^^a„ 

n— \n— n— 771— n+1 



the conditions Ei 



n—O 771— n+1 

A4 in the genus-four configuration just 



• = Eio = determining the endpoints Aq, . 
beyond the secondary caustic have a different form than the relations Ei = ■ ■ ■ = Eiq = (for G = 4) from 
S 13.4.11 In particular, a configuration Ao(a;,i), . . . , A5(x,t) satisfying i?i = • • • = Eiq =■ does not satisfy 
the equations E\ = ■ ■ ■ = Eiq = 0. Nonetheless, these latter conditions also determine implicitly functions 
Xo{x,t), . . . , X5{x, t) satisfying exactly the same Whitham equations. 

Theorem 1. Any continuously differentiable solution Ao(x,t), . . . ,X4{x,t) of Ei = ■ ■ ■ = Eiq = satisfies 
the Whitham equations (|159|l with complex characteristic velocities given by (|160|l in any region of the {x,t)- 
plane where the endpoints are distinct. 

Proof. We complexify the functions En, letting v denote the vector (Ao, . . . , A4, Aq, . . . , X't)'^ of endpoints and 
their complex conjugates, viewed as independent variables. Then, direct calculations show that, by analogy 

with n^ . 



(247) 

and by analogy with H154|l we have 

(248) 
so that 

.dMa 



dvk 2 rj- Vk 



(249) 



dEn 
dvk 



dvk 



R{v)dri 



Ai 



dMo 

dvk 



A! 



Rjri) dy 
V - Vk 



= Jnk{v) 



dvk 



_i 'n-vk 

for n = 1, n — 2, and n = A, where Jnkiy) are exactly the same matrix elements as defined for G = 4 by 
(|155|l . The corresponding calculation for partial derivatives of E^ = V2 is slightly more complicated due to 
the appearance of terms explicitly depending on A2 and A2; however since 



(250) 



fact 



(251) 



d 



2ie°{X2) - 2tt 



drj 



0, 



holds for 71 = 3 as well. Next, by another direct calculation, we find that by analogy with (|156|l . 



dEn 
dvk 



1 ~ 

^En-l 



dEn 



- Vk- 



dvk 
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for 6 < ?? < 10 



so that on a solution we have 

(252) — — = - — = Vf^ — — = Jnk[v)^— , for 5 < ?i < 10 , 

OVk OVk OVk ov 

where again the matrix elements Jnk{v) are precisely the same as defined by (|157f) . Therefore, the elements 
dEn/dvk are assembled into the Jacobian matrix having the factorized form: 

(253) — = J{v) ■ diag 



V \ ovi ovio I 

where the factor 3{v) is exactly the same function of v as in the case studied in S 13.4.11 Noting that the 
partial derivatives dEn/dx and dEn/dt are all identical expressions in terms of v as in S 13.4. II it then follows 
by application of Cramer's rule that 

(254) ^ + c,(z;)^=0, for fc = 1, . . . , 10 , 

where the characteristic speeds are exactly the same functions of v (see (|160ll ) as in § 13.4.11 

To reiterate, the only difference between the derivation here and that in the case when no endpoints lie on 
i^_3 is in the diagonal factor of the Jacobian consisting of the partial derivatives dMo/dvk, which certainly 
are different functions of v than the derivatives dMo/Ovk- However, this diagonal factor cancels out of the 
Cramer's rule formula for the characteristic velocities Ck{v). Note that the condition that the diagonal factor 
is nonsingular is, due to H248() . equivalent to the function R{i])Y(r]) which defines the bands vanishing exactly 
like a square root at the band endpoints. □ 

7. Conclusion 

In this paper we have done the following things. 

• We used a numerical linear algebra method for computing the iV-soliton for several values of as 
large as N = 40, plotted the results and described the phenomena of macrostructure, microstructure, 
and phase transitions (caustic curves) that become more well-defined as A'^ — > oo. 

• We have described a new modification of the asymptotic technique first used in |10| to analyze the 
iV-soliton for large N that gives an improved error estimate of 0{N^^) over the cruder estimate of 
0(Ar-i/3) reported in [TUj . 

• We analyzed the Fourier power spectrum of the A^-soliton in the genus-zero region before the primary 
caustic and related the results to supercontinuum generation of coherent broadband spectra from 
narrow-band sources in optical fibers with weak anomalous dispersion. 

• We computed numerically for the first time the bands of discontinuity of g'{X) in the genus-two 
region beyond the primary caustic and pointed out the phenomenon of a band crossing the locus of 
accumulation of discrete eigenvalues for the A^-soliton. 

• We introduced a new multi-interpolant asymptotic procedure that allows the genus-two ansatz for 
g{X) to be continued through the discrete eigenvalues, proving that this phenomenon is not the cause 
of the secondary caustic. 

• We observed numerically the failure of a new inequality associated with the multi-interpolant ap- 
proach that becomes necessary once a band crosses the discrete spectrum, and we gave evidence that 
this failure is the mathematical mechanism for the secondary caustic phase transition. 

• We have shown that, under the assumption that the new inequality indeed fails as we have observed 
numerically, the secondary caustic represents a change of genus from G = 2 to G = 4 but with new 
topological features of the genus-four ansatz for g{X) that have not been seen before. We also proved 
that despite the unusual appearance of the implicit relations defining the band endpoints in this 
case, the endpoints still satisfy the Whitham equations. 

The nature of the mathematical mechanism for the secondary caustic leads us to hypothesize that there 
could easily be an infinite number of caustic curves in the {x, i)-plane of the A^-soliton in the large N limit. 
The scenario we imagine is the following: as t increases for fixed a; > 0, the new band I2 that we have shown 
opens up upon crossing the secondary caustic drifts toward the imaginary interval [0, iA] of accumulation 
of discrete eigenvalues. When I2 meets this interval, the analysis fails once again until it is rescued (in 
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analogy with the modification we introduced in I3.9|l by adjoining a new contour lobe, -flT-s, enclosing a 
region Z)_5 adjacent to D_3 in which the correct interpolation formula to use is the analogue of (|192|l with 
the term —3id'-'{X) in the exponent replaced by —5i9'^{X). Along K-^ we require a new inequality, namely 
3?((/)(A) — 4i0°(A)) < 0. The region in which this inequality holds eventually pinches off (a tertiary caustic), 
and the situation is resolved by introducing a new band I3 on if- 5 which results in local dynamics described 
by Ricmann thcta functions of genus six. The reader can then imagine that the process repeats again and 
again. The fcth caustic curve is described by eliminating A from the equations 

(255) (j>'{\)-2i{k-l)e"'{\) = 0, di{(f>{X)-2i{k-l)e"{X))^0. 

The formula for (j){X) of course changes each time a caustic is crossed and a new band is added. Thus, a 
harbinger of a developing caustic curve is a recently born band meeting the imaginary axis. This can be 
detected just by observing the wave field as long as (like we see in the numerical constructions of g{X) in 
the genus-two case) the endpoint crossing the imaginary axis has an imaginary part much larger than those 
of all other cndpoints. In this case, the wave field, although locally described by Ricmann theta functions 
of some even genus G > 0, is well-approximated by a complex exponential function g^ikx-u)t)/h (^because the 
genus is "almost" zero), and the real part of the dominant endpoint is proportional to the wavenumber k. 
Thus a sign change of k, which could be detected by local spectral analysis, indicates the dominant endpoint 
crossing the imaginary axis, an event that precedes the formation of a new caustic. 

We want to stress that the mechanism for all of the caustics beyond the primary caustic depends in 
a crucial way on the discrete nature of the spectrum for the iV-soliton, which causes us to reinvent our 
method of analysis each time a band wants to cross the discrete spectrum. Significantly, if we chose to 
begin our analysis by assuming that there is no essential cost in the limit ?i — > for analyzing the formal 
continuum-limit problem for P(A) (in which we simply neglect the difference between the discrete spectral 
measure and its weak limit as — > 00) in place of the exact Riemann-Hilbert problem for P(A), we would 
not need to modify our analysis if a band crosses the spectrum, and we would not need to introduce the new 
inequality 3?(0(A) — 2i9'^{X)) < that leads to the secondary caustic. In other words, the formal continuum- 
limit problem corresponds to a solution ipN{x,t) of the focusing nonlinear Schrodinger equation ^ that 
accurately resembles the iV-soliton ippf{x,t) for {x,t) in the region before the secondary caustic, but not 
beyond. In particular, the solution obtained from the formal continuum-limit problem does not experience 
any phase transition at all when the A^-soliton experiences the secondary caustic. Of course, as the formal 
continuum-limit problem may break the even symmmetry of the A'^-soliton (recall § 13.2.2(1 we have no solid 
grounds for expecting any similarity between 'ipN{x,t) and tpNix^t) for any x and t, and therefore it seems 
to us that wc must therefore regard the correspondence before the secondary caustic as a lucky coincidence. 
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